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Preface 



This volume contains the papers accepted for publication at IPCO X, the Tenth 
International Conference on Integer Programming and Combinatorial Optimiza- 
tion, held in New York City, New York, USA, June 7-11, 2004. The IPCO series 
of conferences presents recent results in theory, computation and applications of 
integer programming and combinatorial optimization. 

These conferences are sponsored by the Mathematical Programming Society, 
and are held in those years in which no International Symposium on Mathemati- 
cal Programming takes place. IPCO VIII was held in Utrecht (The Netherlands) 
and IPCO IX was held in Cambridge (USA). 

A total of 109 abstracts, mostly of very high quality, were submitted. The 
Program Committee accepted 32, in order to meet the goal of having three days 
of talks with no parallel sessions. Thus, many excellent abstracts could not be 
accepted. 

The papers in this volume have not been refereed. It is expected that revised 
versions of the accepted papers will be submitted to standard scientific journals 
for publication. 

The Program Committee thanks all authors of submitted manuscripts for 
their support of IPCO. 



March 2004 



George Nemhauser 
Daniel Bienstock 




Organization 



IPCO X was hosted by the Computational Optimization Research Center 
(CORC) , Columbia University. 

Program Committee 

Egon Balas 
Daniel Bienstock 
Robert E. Bixby 
William Cook 
Gerard Cornuejols 
William Cunningham 
Bert Gerards 
Ravi Kannan 
George Nemhauser, Chair 
William Pulleyblank 
Laurence A. Wolsey 

Organizing Committee 

Daniel Bienstock, Chair 
Garud Iyengar 
Jay Sethuraman 
Cliff Stein 

Sponsoring Institutions 

Bob and Betty Bixby 

IBM 

ILOG 

The Fu Foundation School of Engineering and Applied Science, Columbia Uni- 
versity 

Mathematical Programming Society 




Table of Contents 



Session 1 

Robust Branch-and-Cut-and-Price for the Capacitated Vehicle Routing 

Problem 1 

R. Fukasawa, J. Lysgaard, M. Poggi de Aragao, M. Reis, E. Uchoa, 

R. F. Werneck 

Metric Inequalities and the Network Loading Problem 16 

P. Avella, S. Mattia, A. Sassano 

Valid Inequalities Based on Simple Mixed-Integer Sets 33 

S. Dash, O. Giinliik 

Session 2 

The Price of Anarchy when Costs Are Non-separable and Asymmetric ... 46 

G. Perakis 

Computational Complexity, Fairness, and the Price of Anarchy of the 

Maximum Latency Problem 59 

J.R. Correa, A.S. Schulz, N.E. Stier Moses 

Polynomial Time Algorithm for Determining Optimal Strategies in 

Cyclic Games 74 

D. Lozovanu 

Session 3 

A Robust Optimization Approach to Supply Chain Management 86 

D. Bertsimas, A. Thiele 

Hedging Uncertainty: Approximation Algorithms for Stochastic 
Optimization Problems 101 

R. Ravi, A. Sinha 

Scheduling an Industrial Production Facility 116 

E. Asgeirsson, J. Berry, C.A. Phillips, D.J. Phillips, C. Stein, J. Wein 

Session 4 

Three Min-Max Theorems Concerning Cyclic Orders of Strong Digraphs . 132 

S. Bessy, S. Thomasse 




X 



Table of Contents 



A TDI Description of Restricted 2-Matching Polytopes 139 

G. Pap 

Enumerating Minimal Dicuts and Strongly Connected Subgraphs and 

Related Geometric Problems 152 

E. Boros, K. Elbassioni, V. Gurvich, L. Khachiyan 

Session 5 

Semi-continuous Cuts for Mixed-Integer Programming 163 

I. R. de Farias Jr. 

Combinatorial Benders’ Cuts 178 

G. Codato, M. Fischetti 

A Faster Exact Separation Algorithm for Blossom Inequalities 196 

A.N. Letchford, G. Reinelt, D.O. Theis 

Session 6 

LP-based Approximation Algorithms for Capacitated Facility Location . . 206 
R. Levi, D.B. Shmoys, C. Swamy 

A Multi-exchange Local Search Algorithm for the Capacitated Facility 
Location Problem 219 

J. Zhang, B. Chen, Y. Ye 

Separable Concave Optimization Approximately Equals Piecewise 

Linear Optimization 234 

T.L. Magnanti, D. Stratila 

Session 7 

Three Kinds of Integer Programming Algorithms Based on Barvinok’s 

Rational Functions 244 

J.A. De Loera, D. Haws, R. Hemmecke, P. Huggins, R. Yoshida 

The Path-Packing Structure of Graphs 256 

A. Sebo, L. Szego 

More on a Binary-Encoded Coloring Formulation 271 

J. Lee, F. Margot 

Session 8 

Single Machine Scheduling with Precedence Constraints 283 

J.R. Correa, A.S. Schulz 




Table of Contents 



XI 



The Constrained Minimum Weighted Sum of Job Completion Times 

Problem 298 

A. Levin, G.J. Woeginger 

Session 9 

Near-Optimum Global Routing with Coupling, Delay Bounds, and 

Power Consumption 308 

J. Vygen 

A Flow-Based Method for Improving the Expansion or Conductance of 
Graph Cuts 325 

K. Lang, S. Rao 

All Rational Polytopes Are Transportation Polytopes and All Polytopal 

Integer Sets Are Contingency Tables 338 

J. De Loera, S. Onn 

Session 10 

A Capacity Scaling Algorithm for M-convex Submodular Flow 352 

S. Iwata, S. Moriguchi, K. Murota 

Integer Concave Cocirculations and Honeycombs 368 

A. V. Karzanov 

Minsquare Factors and Maxfix Covers of Graphs 388 

N. Apollonio, A. Sebo 

Session 11 

Low-Dimensional Faces of Random 0/ 1-Poly topes 401 

V. Kaibel 

On Polylredra Related to Even Factors 416 

T. Kirdly, M. Makai 

Optimizing over Semimetric Polytopes 431 

A. Frangioni, A. Lodi, G. Rinaldi 

Author Index 445 




Robust Branch-and-Cut-and-Price for the 
Capacitated Vehicle Routing Problem 



Ricardo Fukasawa 1 , Jens Lysgaard 2 , Marcus Poggi cle Aragao 3 , Marcelo Reis 3 , 
Eduardo Uclroa 4 *, and Renato F. Werneck 5 

1 School of Industrial and Systems Engineering, GeorgiaTech, USA 
rf ukasaw@isye . gatech . edu. 

2 Department of Management Science and Logistics, 

Aarhus School of Business, Denmark 
lys@asb . dk 

3 Departamento de Informatica, PUC Rio de Janeiro, Brazil 
{poggi ,mreis}@inf .puc-rio .br 
4 Departamento de Engenharia de Produgao, 

Universidade Federal Fluminense, Brazil. 
uchoa@producao.uff .br 

5 Department of Computer Science, Princeton University, USA 
rwerneckOcs . princeton . edu 



Abstract. The best exact algorithms for the Capacitated Vehicle Rout- 
ing Problem (CVRP) have been based on either branch-and-cut or La- 
grangean relaxation/column generation. This paper presents an algo- 
rithm that combines both approaches: it works over the intersection 
of two polytopes, one associated with a traditional Lagrangean relax- 
ation over q-routes, the other defined by bound, degree and capacity 
constraints. This is equivalent to a linear program with exponentially 
many variables and constraints that can lead to lower bounds that are 
superior to those given by previous methods. The resulting branch-and- 
cut-and-price algorithm can solve to optimality all instances from the 
literature with up to 135 vertices. This doubles the size of the instances 
that can be consistently solved. 



1 Introduction 

Let G = (V, E) be an undirected graph with vertices V ={0,1,..., n}. Vertex 0 
represents the depot, whereas all others represent clients, each with an associated 
demand d(-). Each edge e £ E has a nonnegative length £(e). Given G and two 
positive integers ( K and C), the Capacitated Vehicle Routing Problem (CVRP) 
consists of finding routes for K vehicles satisfying the following constraints: (i) 
each route starts and ends at the depot, (ii) each client is visited by a single 
vehicle, and (iii) the total demand of all clients in any route is at most G. The 
goal is to minimize the sum of the lengths of all routes. This classical NP-hard 
problem is a natural generalization of the Travelling Salesman Problem (TSP), 

* Corresponding author. 
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and has widespread application itself. The CVRP was first proposed in 1959 by 
Dantzig and Ramser [13] and has received close attention from the optimization 
community since then. 

A landmark exact algorithm for the CVRP, presented in 1981 by Clrristofides, 
Mingozzi and Toth [11], uses a Lagrangean bound from minimum g-route sub- 
problems. A q-route is a walk that starts at the depot, traverses a sequence of 
clients with total demand at most C, and returns to the depot. Some clients may 
be visited more than once, so the set of valid CVRP routes is strictly contained 
in the set of g-routes. The resulting branch-and-bound algorithm could solve 
instances with up to 25 vertices, a respectful size at the time. 

Several other algorithms using Lagrangean relaxation appear in the liter- 
ature. Clrristofides et al. [11] also describe a lower bound based on /c-degree 
center trees, which are minimum spanning trees having degree K < k < 2K 
on the depot, plus 2 K — k least-cost edges. Lagrangean bounds based on K- 
trees (sets of n + K — 1 edges spanning V) having degree 2 K in the depot were 
used by Fisher [14] and by Martinlron, Lucena, and Maculan [24], among others. 
Miller [25] presented an algorithm based on minimum 6-matchings having degree 
2K at the depot and 2 on the remaining vertices. Lagrangean bounds can be 
improved by dualizing capacity inequalities [14,25] and also comb and multistar 
inequalities [24]. 

Another family of exact algorithms stems from the formulation of the CVRP 
as a set partitioning problem by Balinski and Quandt [8] . A column covers a set 
of vertices S with total demand not exceeding C and has the cost of a minimum 
route over {0} U S. Unfortunately, the formulation is not practical because pric- 
ing over the exponential number of columns requires the solution of capacitated 
prize-collecting TSPs, a problem almost as difficult as the CVRP itself. Agarwal, 
Martlrur and Salkin [7] proposed a column generation algorithm on a modified 
set partitioning problem where column costs are given by a linear function over 
the vertices yielding a lower bound on the actual route cost. Columns with the 
modified cost can be priced by solving easy knapsack problems. Hadjiconstanti- 
nou et al. [17] derive lower bounds from heuristic solutions to the dual of the 
set partitioning formulation. The dual solutions are obtained by the so-called 
additive approach, combining the g-route and fc-shortest path relaxations. 

For further information and comparative results on the algorithms mentioned 
above, we refer the reader to the surveys by Toth and Vigo [31,32]. 

Recent research on the CVRP has been concentrated on the polyhedral de- 
scription of the convex hull of the edge incidence vectors that correspond to 
K feasible routes and on the development of effective separation procedures 
[1,3,5,6,12,20,26]. In particular, Araque et al. [4], Augerat et al. [6], Blasum 
and Hochstattler [9], Ralphs et al. [30], Achuthan, Caccetta, and Hill [2] and 
Lysgaard, Letchford, and Eglese [23] describe complete branch-and-cut (BC) al- 
gorithms. These are the best exact methods currently available for the CVRP. 
However, the addition of several elaborate classes of cuts does not guarantee tight 
lower bounds, especially for large values of K (K > 7, say). Closing the resulting 
duality gap usually requires exploring several nodes in the branch-and-cut tree. 
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Even resorting to massive computational power (up to 80 processors running in 
parallel in a recent work by Ralphs [29,30]) several instances with fewer than 80 
vertices, including some proposed more than 30 years ago by Christohdes and 
Eilon [10], can not be solved at all. In fact, branch-and-cut algorithms for the 
CVRP seem to be experiencing a “diminishing returns” stage, where substan- 
tial theoretical and implementation efforts achieve practical results that are only 
marginally better than those of previous works. 

We present a new exact algorithm for the CVRP that seems to break through 
this situation. The main idea is to combine the branch-and-cut approach with the 
g-routes approach (which we interpret as column generation instead of the orig- 
inal Lagrangean relaxation) to derive superior lower bounds. Since the resulting 
formulation has an exponential number of both columns and rows, this leads to 
a branclr-and-cut-and-price (BCP) algorithm. Computational experiments over 
the main instances from the literature show that this algorithm is very consis- 
tent on solving instances with up to 100 vertices. Eighteen open instances were 
solved for the first time. 

The idea of combining column and cut generation to improve lower bounds 
has rarely been used, since new dual variables corresponding to separated cuts 
may have the undesirable effect of changing the structure of the pricing sub- 
problem. However, if cuts are expressed in terms of variables from a suitable 
original formulation, they can be incorporated into the column generation pro- 
cess without disturbing the pricing. We refer to branch-and-bound procedures 
based on such formulations as robust branch- and- cut- and-price algorithms. Poggi 
de Aragao and Uchoa [28] present a detailed discussion on this subject, including 
new reformulation techniques that extend the applicability of robust branch-and- 
cut-and-price algorithms to virtually any combinatorial optimization problem. 
This article on the CVRP is part of a larger effort to demonstrate that these 
methods lead to significant improvements on a wide variety of problems. Ma- 
jor advances have already been reported on two other problems: capacitated 
minimum spanning tree [15] and generalized assignment [27]. 

This article is organized as follows. Section 2 describes the integer program- 
ming formulation we will deal with. Section 3 gives a general description of our 
algorithm, including its two main components: column and cut generation. Fol- 
lowing the work of Irnich and Villeneuve [18] on the CVRP with time windows, 
our column generation procedure eliminates g-routes with small cycles. The sep- 
aration routines are based on the families of inequalities recently discussed by 
Letchford, Eglese, and Lysgaarcl [20,23]. Section 4 presents an empirical analysis 
of our method. Final remarks are made in Sect. 5. 

2 The New Formulation 

A classical formulation for the CVRP [19] represents by Xij the number of times 
a vehicle traverses the edge (i, j) € E. The set of client vertices is denoted by 
V+ = {1, . . . , n}. Given a set S C V + , let d(S) be the sum of the demands of all 
vertices in S, and let S(S) be the cut-set defined by S. Also, let k(S ) = \d,(S)/C~\. 
Define the following polytope in 1R^: 
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Y, x e =2 
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Constraints (1) state that each client is visited once by some vehicle, whereas 
(2) states that K vehicles must leave and enter the depot. Constraints (3) are 
rounded capacity inequalities, which require all subsets to be served by enough 
vehicles. Constraints (4) enforce that each edge not adjacent to the depot is 
traversed at most once (edges adjacent to the depot can be used twice when 
a route serves only one client). The integer vectors x in P\ define all feasible 
solutions for the CVRP. There are exponentially many inequalities of type (3), 
so the lower bound given by 



L\ = min > £ e x e 

xePi ' 
e£E 

must be computed by a cutting plane algorithm. 

Alternatively, a formulation with an exponential number of columns can be 
obtained by defining variables (columns) that correspond to g-routes without 
2-cycles (subpaths i — * j — * i, i yf 0). Restricting the (/-routes to those without 
such cycles improves the formulation and does not change the complexity of the 
pricing [11]. Let Q be an mxp matrix where the columns are the edge incidence 
vectors of all p such (/-routes. Let </J be the coefficient associated with edge e in 
the j-tlr column of Q. Consider the following polytope in JR^, defined as the 
projection of a polytope in ]R P+ ^: 



P2 



Constraints (5) define the coupling between variables x and A. Constraint (6) 
defines the number of vehicles to use. It can be shown that the set of integer vec- 
tors in P 2 also defines all feasible solutions for the CVRP. Due to the exponential 
number of variables A, the lower bound given by 

L 2 = min > £ e x e 
16P2 ' 

e£E 

must be computed using column generation or Lagrangean relaxation. 
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The description of polyhedra associated with column generation or Lagrang- 
ean relaxation in terms of two sets of variables, A and x, used in the definition 
of P 2 , is called Explicit Master in [28]. The main contribution of this article is 
a formulation that amounts to optimizing over the intersection of polytopes P± 
and Pi. The Explicit Master format describes such intersection as follows: 



P 3 = PinP 2 = proj x { 



E Qj ■ E 

3= 1 

P 

E E 

3 = 1 



At 



E x e = 2 
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(2) 
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(4) 
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(5) 
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Ve G E 

V j G p} ■ 



Constraint (6) can be discarded, since it is implied by (2) and (5). Computing 
the improved lower bound 



£3 



min E £ e x 



x£P 3 



e£E 



e 



requires solving a linear program with an exponential number of both variables 
and constraints. A more compact LP is obtained if every occurrence x e in (l)-(4) 
is replaced by its equivalent given by (5). The resulting LP will be referred to 
as the Dantzig-Wolfe Master problem (DWM): 



DWM 
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Capacity inequalities are not the only ones that can appear in the DWM. A 
generic cut E e ^E a e x e > b can be included as Ej=i(E e sB a e9j) ' A j > b. In 
fact, we added all classes of cuts described in [23]. 
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3 The Branch-and-Cut-and-Price Algorithm 



Column Generation. The reduced cost of a column (A variable) is the sum 
of the reduced costs of the edges in the corresponding g-route. Let /i, v, 7 r, 
and u> be the dual variables associated with constraints (8), (9), (10), and (11), 
respectively. The reduced cost c e of an edge e is given by: 



c 



e — 



4 - Mi - Mi - Z TS - LU e 
S\S(S)3e 

4 - v - fij - Z 

S|5(S)3e 



e = {i,j}e£\<5({0}) 
e = {0,j}e<5({0}) . 



An additional generic cut Zj=i(Zee.E a eQj) • 4 — & i n DWM, with dual variable 
a, contributes with the value — a e ■ a to the computation of c e . 

The pricing subproblem consists of finding g-routes (columns) of minimum 
reduced cost. Although this problem is NP-hard in general, it can be solved in 
pseudopolynomial time if all demands are integer [11]. 

The basic data structure is a C x n matrix M. Each entry M (d, v) will 
represent the least costly walk that reaches vertex v using total demand exactly 
d. The entry contains a label consisting of the vertex (v), the cost of the walk, 
and a pointer to a label representing the walk up to the previous vertex. The 
matrix is filled by dynamic programming, row by row, starting with d = 1. For 
each row d , the algorithm goes through each entry v and, for each neighbor 
w of v, evaluates the extension of the walk represented by M(d,v) to w. If 
c(M(d, u)) + C(„ )U ,) < c(M(d + d(w),w)), M(d + d(w),w) is updated. Eventually, 
we will have the most negative walk with accumulated demand at most C that 
arrives at each vertex v. Extending the walk to the depot (whose demand is 
zero), we obtain the corresponding g-route. All negative g-routes thus found 
(there will be at most n) are added to the linear program. There are nC entries 
in the matrix, and each is processed in O(n) time, so the total running time is 
0{n 2 C). 

To strengthen the formulation, we look for s-cycle-free g-routes, for small 
values of s. The algorithm operates as above, using dynamic programming to fill 
a C x n matrix with partial walks. Each entry in the matrix is no longer a single 
label, but a bucket of labels. Bucket M(d,v ) represents not only the cheapest 
s-cycle-free walk with total demand d that ends in v, but also alternative walks 
that ensure that all possible s- vertex extensions from v are considered. Deciding 
which walks to keep is trivial for s = 2 [11], but quite intricate for s > 3. We 
handle this with the approach originally proposed by Irnich and Villeneuve [18] 
(the reader is referred to their paper for details). They show that buckets must 
have size at least s!, so the method quickly becomes impractical. 

We use three heuristics to accelerate the algorithm. They all find only s- 
cycle-free g-routes, but not necessarily the minimum one. Hence, whenever they 
find no negative route, the exact algorithm must verify that none exists. 

The first acceleration technique we use is scaling , which considers g > 1 
units of demand at a time. The exact algorithm is run with a modified demand 
d' v (g) = \d v /g] for every vertex v, and modified capacity C' = \C/g\ for the 
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vehicles. The running time will be proportional to C instead of C , but the 
solutions found will still be valid in the original setting. 

The second technique is sparsification. Intuitively, small edges in the original 
graph are more likely to appear in the solution. With that in mind, we pre- 
compute a set of five disjoint spanning forests using an extension of Kruskal’s 
algorithm similar to the “edge-greedy heuristic” described in [34] . By considering 
only edges in those trees, the dynamic programming will be limited to a restricted 
set of neighbors when processing each vertex. 

The third acceleration is bucket pruning , which consists of storing only s 
(instead of s!) labels per bucket. To minimize the number of valid s-extensions 
blocked, we keep only labels that differ significantly from each other. 

The algorithm starts using all three acceleration techniques. Whenever it 
fails to find a negative g-route, it turns one of them off. If the algorithm fails 
with no acceleration, we can safely declare that no negative g-route exists. On 
the other hand, if only a subset of the heuristics is in use and they do find a 
negative g-route, other heuristics are turned back on. 

Cut Generation. At first, the LP formulation includes only degree constraints; 
violated cuts are added as needed. Besides the well-known rounded capacity cuts 
(10) and bound cuts (11), we also use framed capacity, strengthened comb, mul- 
tistar, partial multistar, generalized multistar and lrypotour cuts, all presented 
in [23]. We use the CVRPSEP package [22] to separate those cuts. Since they 
are defined in terms of the x* variables, not readily available in the DWM for- 
mulation, we convert A* variables into the corresponding x* variables as needed. 

Branching Rule. The algorithm starts by adding columns to the formulation 
until there are none with negative reduced costs, then it looks for all violated 
cuts. This is repeated until both column generation and separation fail. If the 
node cannot be fathomed at this point, we branch. 

The branching rule is an adaptation of the rule used in [23]. We choose as 
branching candidates sets S such that 2 < x*(S(S )) < 4 and branch by imposing 
the disjunction (x(S(S)) = 2) V (a;(i5(S')) > 4). 

We use strong branching to select which set to branch on. Since computing 
the actual lower bound for each child node can be very time-consuming, we 
estimate it by performing a small number of column generation iterations. Only 
p candidate sets (5 < p < max{10 — depth, 5}) are evaluated. This limits the 
time spent on strong branching, especially when the depth is high. Priority is 
given to sets with smaller values of |x*(<5(5)) — 2.7|/d(S). We use 2.7 because 
constraint x(S(S)) = 2 tends to have a greater impact on the LP relaxation 
than x(S(S)) > 4, leading to an imbalance in the branch-and-bound tree. Values 
closer to 2 than 4 increase the impact of imposing x(8(S)) > 4. 

Node Selection and Primal Bounds. We used as initial primal bound the 
best previously known solution plus one. The node selection strategy chosen was 
depth first search, because it requires less memory. 
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Dynamic Choice of Column Generation. It usually pays off to combine 
the strengths of cut and column generation into a single algorithm. However, in 
some instances, the increase in the lower bound is not worth the additional time 
spent. In these cases pure branch-and-cut performs better than branch-and-cut- 
and-price. We modified our algorithm to handle such cases. In the root node, we 
start as in a pure branch-and-cut, without any column generation. After that, 
the root node is solved again with a branch-and-cut-and-price using the column 
generation with s = 2. Then we solve it once again, changing s to 3. For the rest 
of the branch-and-bound tree, the algorithm picks the strategy with the best 
balance between running time and bound quality. 

4 Computational Experiments 

We tested our algorithm on most instances available at www . branchandcut . org. 
All instances with up to 135 vertices were solved to optimality, eighteen of them 
for the first time. The tests were run on a Pentium IV, 2.4 GHz, with 1GB of 
RAM. Throughout this section, we will refer to the branch-and-cut-and-price al- 
gorithm with dynamic choice of column generation as Dyn-BCP and the branch- 
and-cut-and-price with pre-determined column generation as BCP. 

Tables 1 and 2 detail the results obtained by Dyn-BCP. Columns Root LB 
and Root Time give the lower bound and total CPU time of the root node 
of the branch-and-bound tree. The running time includes the dynamic choice of 
column generation and strong branching. The next column represents the size 
s of the cycles eliminated by the column generation procedure indicates 
that column generation was not used). Tree Size represents the total number 
of branch-and-bound nodes explored, and Total Time is the total CPU time 
spent by the algorithm. Finally, Prev UB and Our UB indicate the best 
solution value available in the literature and the one obtained by our procedure. 
Values in bold indicate proven optimality. 

Table 3 compares several different lower bounds obtained for the CVRP. 
Lower bounds LI, L2 and L3 (defined in Sect. 2) are presented in the first three 
columns. Column LLE03 shows the lower bounds obtained by [23], which are the 
best available in the literature. The next three columns contain the lower bounds 
obtained at the root node of the BCP with s-cycle elimination (for s = 2,3,4). 
Column OPT is the optimal solution value and T(3) is the total CPU time (in 
seconds) spent at the root node of BCP for s = 3, the value most commonly 
used. On average, the algorithm for s = 4 was 61% slower than for s = 3, which 
in turn was 12% slower than for s = 2. The last line shows the average gap of 
each bound, considering only instances with a known optimum. 

Table 3 shows that major improvements can be obtained by using the BCP 
approach instead of pure branch-and-cut. Also, eliminating longer cycles often 
increases the lower bounds significantly. 

The greatest improvements usually happen when the ratio n/K is smallest. In 
such cases, columns without small cycles are very close to being actual routes, so 
one can expect the duality gap to be reduced. Figure 1, created from the instances 
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Table 1 . Results of the Dyn-BCP algorithm for the A and B instances. 



Instance 


Root 

LB 


Root 
Time (s) 


s 


Tree 

Size 


Total 
Time (s) 


Prev. 

UB 


Our 

UB 


A-n37-k5 


664.8 


16 


- 


8 


19 


669 


669 


A-n37-k6 


932.6 


30 


3 


74 


379 


949 


949 


A-n38-k5 


716.5 


12 


- 


52 


26 


730 


730 


A-n39-k5 


816.6 


107 


3 


11 


167 


822 


822 


A-n39-k6 


822.8 


39 


3 


11 


98 


831 


831 


A-n44-k6 


934.8 


52 


2 


6 


90 


937 


937 


A-n45-k6 


938.1 


52 


3 


11 


170 


944 


944 


A-n45-k7 


1139.3 


88 


3 


26 


331 


1146 


1146 


A-n46-k7 


914.0 


63 


2 


3 


92 


914 


914 


A-n48-k7 


1069.1 


72 


3 


8 


166 


1073 


1073 


A-n53-k7 


1003.9 


138 


3 


16 


611 


1010 


1010 


A-n54-k7 


1153.9 


125 


3 


90 


1409 


1167 


1167 


A-n55-k9 


1067.4 


32 


3 


7 


84 


1073 


1073 


A-n60-k9 


1344.4 


161 


3 


224 


3080 


1354 


1354 


A-n61-k9 


1022.5 


108 


3 


121 


1883 


1034 


1034 


A-n62-k8 


1280.4 


722 


3 


101 


3102 


1290 


1288 


A-n63-k9 


1607.0 


238 


3 


49 


1046 


1616 


1616 


A-n63-kl0 


1299.1 


136 


3 


387 


4988 


1315 


1314 


A-n64-k9 


1385.3 


265 


3 


648 


11254 


1402 


1401 


A-n65-k9 


1166.5 


154 


3 


17 


516 


1174 


1174 


A-n69-k9 


1141.4 


289 


3 


391 


7171 


1159 


1159 


A-n80-kl0 


1754.0 


1120 


3 


153 


6464 


1763 


1763 


B-n38-k6 


800.2 


10 


- 


14 


14 


805 


805 


B-n39-k5 


549.0 


3 


- 


1 


3 


549 


549 


B-n41-k6 


826.4 


13 


- 


8 


18 


829 


829 


B-n43-k6 


733.7 


13 


- 


74 


29 


742 


742 


B-n44-k7 


909.0 


9 


- 


1 


9 


909 


909 


B-n45-k5 


747.5 


10 


- 


19 


16 


751 


751 


B-n45-k6 


677.5 


224 


3 


3 


279 


678 


678 


B-n50-k7 


741.0 


5 


- 


3 


6 


741 


741 


B-n50-k8 


1295.0 


97 


3 


287 


2845 


1312 


1312 


B-n51-k7 


1025.2 


16 


- 


83 


46 


1032 


1032 


B-n52-k7 


745.8 


7 


- 


9 


9 


747 


747 


B-n56-k7 


704.0 


15 


- 


9 


22 


707 


707 


B-n57-k7 


1149.2 


76 


- 


133 


168 


1153 


1153 


B-n57-k9 


1596.0 


61 


3 


15 


193 


1598 


1598 


B-n63-kl0 


1479.4 


231 


- 


501 


682 


1496 


1496 


B-n64-k9 


859.3 


70 


- 


7 


86 


861 


861 


B-n66-k9 


1307.5 


145 


3 


126 


1778 


1316 


1316 


B-n67-kl0 


1024.4 


218 


- 


327 


568 


1032 


1032 


B-n68-k9 


1263.0 


260 


3 


5912 


87436 


1275* 


1272 


B-n78-kl0 


1215.2 


193 


3 


90 


1053 


1221 


1221 



* An earlier version of our algorithm [16] found a solution of 1272 but could not prove 
its optimality. Using that information, K. Wenger [33] solved this instance. 
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Table 2. Results of the Dyn-BCP algorithm for the E, F, M and P instances. 



Instance 


Root 

LB 


Root 
Time (s) 


s 


Tree 

Size 


Total 
Time (s) 


Prev. 

UB 


Our 

UB 


E-nl3-k4 


247.0 


0 


- 


1 


0 


247 


247 


E-n22-k4 


375.0 


0 


- 


1 


0 


375 


375 


E-n23-k3 


569.0 


0 


- 


1 


0 


569 


569 


E-n30-k3 


533.3 


7 


- 


6 


7 


534 


534 


E-n31-k7 


378.5 


4 


2 


2 


6 


379 


379 


E-n33-k4 


834.5 


14 


- 


5 


15 


835 


835 


E-n51-k5 


518.2 


51 


- 


8 


65 


521 


521 


E-n76-k7 


670.0 


264 


2 


1712 


46520 


682 


682 


E-n76-k8 


726.5 


277 


2 


1031 


22891 


735 


735 


E-n76-kl0 


817.4 


354 


3 


4292 


80722 


830 


830 


E-n76-kl4 


1006.5 


224 


3 


6678 


48637 


1021 


1021 


E-nl01-k8 


805.2 


1068 


3 


11622 


801963 


815 


815 


E-nl01-kl4 


1053.8 


658 


3 


5848 


116284 


1071 


1067 


F-n45-k4 


724.0 


8 


- 


1 


8 


724 


724 


F-n72-k4 


236.4 


70 


- 


42 


121 


237 


237 


F-nl35-k7 


1159.9 


6618 


- 


25 


7065 


1162 


1162 


M-nlOl-klO 


820.0 


119 


- 


1 


119 


820 


820 


M-nl21-k7 


1031.1 


5594 


3 


40 


25678 


1034 


1034 


P-nl6-k8 


449.0 


1 


2 


3 


1 


450 


450 


P-nl9-k2 


212.0 


1 


- 


1 


1 


212 


212 


P-n20-k2 


213.0 


1 


- 


9 


1 


216 


216 


P-n21-k2 


211.0 


1 


- 


1 


1 


211 


211 


P-n22-k2 


216.0 


2 


- 


2 


2 


216 


216 


P-n22-k8 


603.0 


3 


2 


1 


3 


603 


603 


P-n23-k8 


529.0 


18 


2 


1 


18 


529 


529 


P-n40-k5 


456.9 


28 


- 


5 


34 


458 


458 


P-n45-k5 


506.6 


59 


3 


11 


194 


510 


510 


P-n50-k7 


551.5 


79 


3 


7 


143 


554 


554 


P-n50-k8 


616.3 


102 


3 


1084 


9272 


649 


631 


P-n50-kl0 


689.3 


50 


3 


65 


304 


696 


696 


P-n51-kl0 


735.2 


35 


3 


22 


105 


741 


741 


P-n55-k7 


557.9 


90 


2 


450 


4649 


568 


568 


P-n55-k8 


579.8 


42 


2 


196 


1822 


588 


588 


P-n55-kl0 


681.4 


107 


3 


1556 


9076 


699 


694 


P-n55-kl5 


972.8 


251 


3 


398 


1944 


993 


989 


P-n60-kl0 


738.9 


126 


3 


52 


570 


756 


744 


P-n60-kl5 


962.8 


118 


3 


76 


442 


1033 


968 


P-n65-kl0 


786.0 


159 


3 


23 


422 


792 


792 


P-n70-kl0 


814.5 


292 


3 


1752 


24039 


834 


827 


P-n76-k4 


588.8 


363 


- 


59 


572 


593 


593 


P-n76-k5 


616.8 


273 


- 


3399 


14546 


627 


627 


P-nl01-k4 


678.5 


1055 


- 


23 


1253 


681 


681 
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Table 3. Comparison of lower bounds. 



Instance 


LI 


L2 


Lower Bounds 
L3 LLE03 s = 2 


s = 3 


s = 4 


OPT T(3) 


A-n53-k7 


996.6 


978.5 


1002.2 


998.7 


1001.9 


1003.8 


1004.1 


1010 


29 


A-n54-k7 


1130.7 


1114.0 


1150.0 


1135.3 


1150.9 


1153.6 


1156.1 


1167 


26 


A-n55-k9 


1055.9 


1025.4 


1066.4 


1058.3 


1066.4 


1067.3 


1067.6 


1073 


12 


A-n60-k9 


1316.5 


1305.6 


1341.6 


1319.6 


1341.6 


1344.4 


1345.3 


1354 


27 


A-n61-k9 


1004.8 


996.8 


1018.6 


1010.2 


1018.9 


1022.7 


1023.3 


1034 


19 


A-n62-k8 


1244.1 


1222.7 


1274.1 


1251.7 


1273.7 


1280.5 


1281.1 


1288 


36 


A-n63-k9 


1572.2 


1564.8 


1603.5 


1580.7 


1603.6 


1607.0 


1608.9 


1616 


34 


A-n63-kl0 


1262.2 


1267.4 


1294.5 


1266.6 


1294.4 


1299.1 


1302.5 


1314 


28 


A-n64-k9 


1340.1 


1353.3 


1378.9 


1351.6 


1379.0 


1385.3 


1389.2 


1401 


39 


A-n65-k9 


1151.1 


1133.0 


1163.4 


1155.2 


1164.4 


1166.5 


1168.3 


1174 


31 


A-n69-k9 


1108.9 


1113.2 


1138.4 


1114.4 


1139.3 


1141.4 


1144.1 


1159 


87 


A-n80-kl0 


1699.9 


1712.2 


1749.7 


1709.6 


1749.6 


1753.9 


1755.6 


1763 


940 


B-n50-k7 


740.0 


664.8 


741.0 


741.0 


741.0 


741.0 


741.0 


741 


8 


B-n50-k8 


1279.2 


1217.5 


1291.8 


1281.1 


1291.8 


1295.2 


1302.0 


1312 


20 


B-n51-k7 


1024.6 


918.4 


1025.9 


1025.6 


1026.0 


1026.4 


1026.6 


1032 


21 


B-n52-k7 


745.0 


640.2 


746.3 


746.0 


746.4 


746.4 


746.4 


747 


17 


B-n56-k7 


703.4 


606.9 


704.5 


705.0 


704.5 


705.0 


705.0 


707 


19 


B-n57-k7 


1148.6 


1058.5 


1150.9 


1150.1 


1151.1 


1151.6 


1152.3 


1153 


126 


B-n57-k9 


1586.7 


1511.5 


1595.2 


1589.2 


1595.2 


1595.8 


1596.1 


1598 


16 


B-n63-kl0 


1478.9 


1418.4 


1484.2 


1481.0 


1484.3 


1487.2 


1487.3 


1496 


32 


B-n64-k9 


858.5 


769.3 


860.1 


860.5 


860.9 


861.0 


860.5 


861 


44 


B-n66-k9 


1295.2 


1223.1 


1302.6 


1298.5 


1303.7 


1307.5 


1308.5 


1316 


58 


B-n67-kl0 


1023.8 


984.5 


1026.4 


1024.8 


1026.4 


1026.9 


1027.2 


1032 


25 


B-n68-k9 


1256.8 


1163.9 


1261.5 


1258.1 


1261.6 


1262.9 


1263.5 


1272 


42 


B-n78-kl0 


1202.3 


1124.5 


1212.5 


1205.6 


1212.7 


1214.8 


1215.7 


1221 


63 


E-n51-k5 


514.5 


512.9 


518.0 


519.0 


519.1 


519.1 


519.4 


521 


23 


E-n76-k7 


661.4 


663.3 


668.4 


666.4 


670.3 


670.8 


670.9 


682 


73 


E-n76-k8 


711.2 


716.7 


725.1 


717.9 


726.5 


726.8 


727.0 


735 


72 


E-n76-kl0 


789.5 


811.4 


816.5 


799.9 


817.3 


817.4 


817.7 


830 


36 


E-n76-kl4 


948.1 


999.6 


1004.8 


969.6 


1004.9 


1006.5 


1007.5 


1021 


15 


E-nl01-k8 


796.4 


786.4 


801.8 


802.6 


804.6 


805.2 


805.2 


815 


320 


E-nl01-kl4 


1008.3 


1045.1 


1051.6 


1026.9 


1052.2 


1053.9 


1054.3 


1067 


69 


M-nlOl-klO 


819.5 


798.1 


820.0 


820.0 


820.0 


820.0 


820.0 


820 


31 


M-nl21-k7 


1009.7 


1013.0 


1030.9 


1017.4 


1031.2 


1031.8 


1032.0 


1034 


1257 


P-n50-k8 


596.9 


612.5 


615.3 


602.1 


615.8 


616.3 


617.1 


631 


13 


P-n55-kl0 


646.7 


677.2 


680.1 


662.1 


680.2 


681.4 


681.4 


694 


6 


P-n55-kl5 


895.1 


963.0 


967.5 


906.7 


967.7 


973.1 


973.2 


989 


7 


P-n60-kl0 


708.3 


733.5 


737.2 


718.4 


737.7 


739.0 


739.3 


744 


8 


P-n60-kl5 


903.3 


955.6 


961.2 


929.8 


961.7 


962.9 


963.5 


968 


4 


P-n65-kl0 


756.5 


779.8 


785.2 


767.1 


785.6 


786.0 


787.1 


792 


13 


P-n70-kl0 


786.9 


808.3 


812.7 


795.6 


813.6 


814.6 


814.9 


827 


23 


AvgGAP 


2.92% 


4.76% 


1.02% 


2.26% 


0.97% 


0.82% 


0.74% 


— 


— 
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n/K 



Fig. 1 . Gap improvement of BCP with s = 3 (s3) with respect to [23] (lie). The 
graph represents l [f e as a function of j/. Circles represent instances solved to 
optimality for the first time and crosses the ones already closed. 



of Table 3 with a known optimal solution, illustrates this observation. The hori- 
zontal axis represents the ratio n/K. The value plotted is (s3 — lie)/ (opt — lie), 
where s3 is the bound of the BCP with s = 3, lie is the LLE03 bound and opt 
is the optimal solution value. This ratio represents the reduction of the dual- 
ity gap relative to the LLE03 gap. The value is zero whenever both methods 
produce the same result; positive values mean that BCP gaps are smaller, and 
negative values favor LLE03. Circles represent instances that were solved by our 
algorithm for the first time, while crosses indicate instances that were already 
closed. Only one instance (B-n56-k7) had a negative value (—0.009). The gap 
was often reduced by more than 50%, enough to close several open instances. 

5 Conclusion 

This paper has shown that a careful combination of column and cut generation 
leads to a very effective algorithm for the CVRP. We consider the three main as- 
pects that can lead to further improvements: cut generation, column generation, 
and the integration between them. 

Consider cutting planes. Rounded capacity cuts are highly effective in our 
BCP — the efficient heuristic separation for them contributes significantly to the 
good overall performance. We did try all other families of cuts for the CVRP 





Robust Branch-and-Cut-and-Price for the CVRP 



13 



known to be effective in a pure BC: framed capacities, generalized capacities, 
strengthened combs, multistars, and extended hypotours. Surprisingly, however, 
their effect was quite modest in practice — only on E-nl01-k8 the root bound 
increased by more than two units. Any improvement to the cutting part of our 
BCP will therefore require new families of inequalities. Unfortunately, looking 
for facet-defining inequalities may not be enough within the BCP framework, 
since a facet of the CVRP polyhedron can also be a facet of the polyhedron 
induced by column generation over the g-routes. 6 A new family of inequalities 
will ouly be effective in our context if it significantly cuts polyhedron P 2 . 

In the short term, improvements are more likely to come from column gener- 
ation. The goal is to devise efficient ways of pricing over a set of columns as close 
to cycle-free g-routes as possible. The s-cycle elimination approach we chose was 
reasonably effective, at least for s < 4. But other alternatives to further restrict 
the g-routes can also be tried. One could, for instance, forbid g-routes that visit 
vertices in a small set S more than once. This set could be dynamically chosen 
in order to eliminate g-routes with long cycles that are particularly harmful to 
bound quality. 

In principle, several of the combinatorial relaxations of the CVRP — such as 
fc-degree center trees, AT-trees, or 6-matchings — could also be used instead of (or 
even in addition to) g-routes. We think they are unlikely to have a noticeable 
impact, for two reasons. First, among the relaxations mentioned above, only 
g-routes take client demands and vehicle capacities into account. These “nu- 
merical” elements are precisely what makes the CVRP much harder in practice 
than a “pure” graph problem such as the TSP. It seems that known families of 
inequalities for the CVRP, mostly inspired by previous work on the TSP, cannot 
cope well with this numerical aspect. The inequalities implicitly given by P 2 can 
do better. The second reason for not using the alternatives above is that only g- 
routes lead to a pricing problem that is superpolynomial, but still practical. The 
polylreclra associated with polynomial pricing problems can be fully described 
and efficiently separated, which usually makes a pure BC (instead of BCP) a 
faster alternative. For instance, we can separate over the 6-matching polyhedron 
in a very efficient way [21]. Separation over the g-route polyhedron, on the other 
hand, is possible only implicitly, with column generation. 

Even if we leave column and cut generation as they are, we can still accelerate 
the algorithm through better integration. The task of designing, coding and 
fine-tuning a high performance branch-and-cut or branclr-and-price procedure is 
complex, requiring experience to choose among several possible strategies, and, 
often, some computational tricks. BCP multiplies the possibilities and, therefore, 
the difficulties. We believe that there is still plenty of room for such improvements 
in our BCP. 



A. Letchford proved that generalized large multistar inequalities do not cut the 
polyhedron induced by cycle- free g-routes (personal communication). Although this 
is not necessarily true when using g-routes without s-cycles, the heuristic separation 
in our BCP never found a violated GL multistar. 
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Abstract. Given a simple graph G(V,E) and a set of traffic demands 
between the nodes of G, the Network Loading Problem consists of in- 
stalling minimum cost integer capacities on the edges of G allowing rout- 
ing of the traffic demands. 

In this paper we study the Capacity Formulation of the Network Load- 
ing Problem, introducing the new class of the Tight Metric Inequalities, 
that completely characterize the convex hull of the integer feasible solu- 
tions of the problem. We present separation algorithms for Tight Metric 
Inequalities and a cutting plane algorithm, reporting on computational 
experience. 



1 Introduction 

Let G(V, E) be a graph without loops and multiple edges. A set of traffic demands 
have to be routed between the nodes of G. Let dij denote the demand between 
nodes i and j and let D( V, H) be the (directed) demand graph, where H = {in £ 
V x V : d^ > 0} and d : H -► R+. 

A capacity must be installed on the edges of G in integer multiples of a base 
unit. A set of capacities X = {diij : ij € E] is said feasible if all the demands 
can be routed simultaneously with the constraint that the total flow traversing 
ij £ E does not exceed x-ij . 

Let Cij be the cost of installing a unit capacity on the edge ij £ E and 
let c(X) = J^ijeE denote the cost of X. The Network Loading Problem, 
denoted by NL(c, G , d), is the problem of determining a set of feasible capacities 
X minimizing c(X). Let p(c,G,d) denote the value of the optimal solution of 
NL{c, G , d). 

We refer to the Asymmetric Network Loading Problem if a directed graph 
G(V, A) is used. 



D. Bienstock and G. Nemhauser (Eds.): IPCO 2004, LNCS 3064, pp. 16—32, 2004. 
(c) Springer- Verlag Berlin Heidelberg 2004 
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Let O = {o € V : Y doi > 0} be the set of the source nodes of G. By 
iev 

associating a commodity with each source node o£ O , we can write the Flow 
Formulation , where /° denotes the amount of commodity o traversing the edge 
ij: 



min CijXij 




ijeE 




E Pij~ E % = -<%> o £ 0,i £ V 


(1) 


j.ijeE j.ijeE 




E(tf 


(2) 


oGO 




fij > 0, o £ O, ij £ E 


(3) 


Xij £ Z + , ij £ E 


(4) 



Flow constraints (1) guarantee that the traffic demand d° is routed from the 
source node o to each node i £ V. Capacity constraints (2) ensure that the total 
flow traversing ij does not exceed the installed capacity X{j . 

By projecting flow variables, the Network Loading Problem can be formu- 
lated in the space defined by the capacity variables x ij. 

Let i, j £ V and let 77^- denote the family of the paths between i and j on 

I E 1 1 

G. The function p : E — > IR 1 ^ defines a metric on G if Y dhk ~ dij > 0, for 

hk^Pij 

each Pij £ Ilij and ij £ E [ 31 ], and we refer to the set Met(G) = {/a £ : 

Y dhk — l^ij > 0 ,Pij £ IIij,ij £ £} as the metric cone. Let 1} denote the 

hk£Pij 

shortest distance between nodes i and j according to the weights p. 

By using the feasibility condition for multicommodity flows introduced in 
[30] [34] [28], the Capacity Formidation for the Network Loading Problem is: 

min E Cij %ij 

ijeE 

px > Id, p. £ Met(G) (5) 

x e zlf 1 

where inequalities (5) are referred to as Metric Inequalities. 

Let px > Id be any Metric Inequality. The Capacity Formulation can be 
strengthened by studying the Knapsack Polyhedron KN AP(p) = Conv {x £ 

I FI 

ZV : fix > Id}. Valid inequalities for KNAP(p) can be derived by multiplying 
both the sides by A and then rounding-up the l.lr.s. coefficients and the r.lr.s.. 

Let p £ Z^L Inequalities of the form px > \td ], are referred to as Rounded 
Metric Inequalities. Let S, T be two disjoint subsets of V defining a partition and 
let (S : T) = {ij : i £ S,j £ T). Special kinds of Rounded Metric Inequalities 
are the Cut Inequalities x(S : T) > |"d(5 : T) + d{T : 5)]. 
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The Network Loading Problem has received a considerable amount of atten- 
tion in the literature [3]. 

Magnanti, Mirchandani and Vachani [32] [33] studied basic polyhedral prop- 
erties of the problem and introduced three basic classes of valid inequalities: 
Cut Inequalities, 3-Partition Inequalities and Arc Residual Capacity Inequali- 
ties. Atamturk and Raj an [2] presented a linear-time separation algorithm for 
Arc Residual Capacity Inequalities. 

Barahona [4] presented a cutting plane algorithm based on Cut Inequali- 
ties, to provide lower bounds for the Network Loading Problem. The separation 
problem for Cut Inequalities is formulated as a Max-Cut problem. Barahona 
also introduced Spanning Tree Inequalities to enforce connectivity of feasible so- 
lutions. The solution of the relaxed formulation is used to set lower bounds for 
the capacities and then a Branch-and-Bound based on the Flow Formulation is 
run to provide an upper bound. 

Dahl and Stoer [21] studied a Network Loading Problem with survivability 
requirements. They present a cutting plane algorithm where violated Metric 
Inequalities are generated by solving a LP by a column generation approach. 

Bienstock and Giinliik [12] investigated the polyhedral structure of the Flow 
Formulation for the asymmetric version of the problem and developed a cutting 
plane algorithm. They derived the Flow Cutset Inequalities as a class of Mixed- 
Integer Rounding Inequalities. Chopra, Gilboa and Sastry [16] generalized the 
Flow Cutset Inequalities introduced in [12]. 

Bienstock, Chopra, Giinliik and Tsai [10] compared cutting plane algorithms 
based, respectively, on the Flow and on the Capacity Formulation. Giinliik [25] 
introduced new families of Mixed-Integer Rounding Inequalities. 

Atamturk [1] studied the polyhedral properties of the Flow Formulation for 
the special Network Loading Problem over a ’’cutset”, defined, given a nonempty 
partition {U\, U 2 ) of V , as the set of arcs directed from U\ to U 2 and vice-versa. 

Solving the LP-relaxation of the Flow Formulation for large instances, re- 
quires a huge amount of memory. Crainic, Frangioni and Gendron [17] [18] solved 
the lagrangian relaxation of the Flow Formulation by the Bundle Method. Bien- 
stock [6] [7] [8] has developed a potential reduction approach to solve large-scale 
block structured linear programs. The algorithm turned out to be very effective 
to tackle large-scale Network Loading instances. 

Holmberg and Yuan [27] proposed a Branch-and-Bound based on a lagrangian 
relaxation of the Flow Formulation for the case where Xij £ {0,1}, ij £ E. Local 
search heuristics have been investigated in [20] [19]. 

In this paper we study the Capacity Formulation of NL(c, G, d), introducing 
the new class of the Tight Metric Inequalities, that completely characterize the 
convex hull of the integer feasible solutions of the problem. We present separation 
algorithms for Tight Metric Inequalities and a cutting plane algorithm, reporting 
on computational experience. 

The remainder of the paper is organized as follows. In section 2 we intro- 
duce the Tight Metric Inequalities. In section 3 we address several separation 
algorithms for Metric and Tight Metric Inequalities. In section 4 we outline 
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the cutting plane algorithm. Finally in section 5 we report on computational 
experience with two sets of test instances. 

2 Tight Metric Inequalities 

Let P{G,d ) = conv {x € : x is feasible} be the Network Loading Polyhe- 

dron. In this section we introduce the family of the Tight Metric Inequalities 
and we prove that they completely characterize P(G, d). 

We first show that the left-hand-side of any non-redundant valid inequality 
defines a metric. 

Theorem 1 . Let ax > b be any valid inequality for P(G,d). There exists /i £ 
Met{G) such that: 

i) fix > b is valid 
n ) fXij _ o. j j . ij £ E 

Proof. Let i and j be two distinct nodes of V and let a(P™ m ) denote the length 
of the shortest path between i and j according to the weights a. Let pij = 
a{P™ q,n ) < atj , ij £ E. It follows that p, £ Met(G ) and the inequality fix > b 
dominates ax > b. 

Suppose that fix > b is not valid and let i be a feasible solution with the 
property that f ix < b. Without loss of generality, we can assume that there is 
only one edge hk with phk = o-(P]fjf n ) < ahk- 

We can build a new feasible solution x , defined as follows: 

(x i:i ij ^ hk, ij i Ph k j n 

Xij = < Xij + x hk ij £ PfH m , 

[ 0 ij = hk 

with the property that ax = \xx < b. It follows that x is feasible and ax > b is 
violated, contradiction. 

Let /.i £ Met(G) and let p(/x, G, d) be the optimal solution of the Network 
Loading Problem where the vector /i replaces c in the objective function, i.e. 
p(/j.,G,d) = min {/./a: : x £ 1 } + ,x is feasible}. Any valid inequality fix > b for 
P(G, d) is dominated by fix > p(p 1 G, d). We refer to any inequality of the form 
px > p(p, G, d) as a Tight Metric Inequality. 

Example 1 . Let G(V,E ) be a complete graph on 3 nodes, V = { 1 , 2 , 3 }, and 
let d\2 = di3 = ^23 = 1 . 2 . From the integer metric p±2 = P13 = P23 = 1 
we get the Metric Inequality X12 + £13 + X23 > 3 . 6 , that can be strengthened 
to the Rounded Metric Inequality £12 + £13 + £23 > 4 . But it can be easily 
checked that p(p, G,d ) =5 and the corresponding Tight Metric Inequality is 
£12 + £13 + £23 > 5 . 

The following result emphasizes the crucial role of the Tight Metric Inequal- 
ities. 
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Corollary 1. The Tight Metric Inequalities px > p(p,G,d), fi € Met(G), com- 
pletely describe P(G,d). 

A metric p £ Met(G) is extreme if p is an extreme ray of the metric cone 
Met(G), i.e. if p cannot be derived as a conic combination of two other distinct 
metrics. The following property holds if p £ Zl s l is extreme. 

Theorem 2. If p £ Zl E l is an extreme ray of Met(G), then p(p,G,d) = \£d\. 

Proof. Let A be a m x n matrix and let Pi = conv{x £ Z n : Ax > b}. 

Let u £ R™ be a vector of multipliers. Starting from the linear system Ax > b , 
the Chvatal-Gomory procedure produces a valid inequality for Pi of the type 
puA]a; > \ub] and adds it to Ax > b. It is a well known result [36] that any valid 
inequality for Pi can be derived by applying the Chvatal-Gomory procedure a 
finite number of times. 

Now let us turn to our problem, where Ax > b is a set of Metric Inequalities. 
Let 7 = \uA \ , 5 = [ ub] and let jx > 6 be the inequality derived by the Chvatal- 
Gomory procedure. It can be easily proved that 7 € Mef(G). It follows that at 
any stage of the procedure the system A includes only inequalities of the form 
72: > b , with 7 £ Met(G). 

Let p an extreme ray of A Iet(G) and let px > k be valid for P(G, d). Since p 
cannot be obtained as a conic combination of two or more other distinct metrics, 
it follows that px > k can be derived applying the Chvatal-Gomory procedure 
with a vector of multipliers u where the multiplier associated with the px > Id 
itself is 1 and the multipliers associated with the other inequalities are null and 
then k = \Id\. 

Let S and T be a partition of V and let (S : T) = {ij £ E : i £ S, j £ T} 
be a cut on G. Deza and Laurent [22] proved that the metrics defined by cuts, 
i.e. pij = 1 for ij £ (S : T) and ptj = 0 for ij £ A\{S : T) are the only 
extreme metrics with p £ Zl E L It follows that Cut Inequalities x(S : T) > \d(S : 
T) + d(T : 5)] are Tight Metric Inequalities. 

3 Separation Algorithms 

Separation of Tight Metric Inequalities is a difficult task, as even computing 
the r.lr.s p(p,G,d) for a given inequality px > b is NP-harcl. We address the 
separation problem by a two-stage heuristic: i) we look for a violated inequality 
ax > b and then ii) we turn it into the px > 6, and apply a shrinking procedure to 
compute p(p, G, d) to get a Tight Metric Inequality of the form px > p(p , G, d). 

The remainder of this section is organized as follows. In subsection 3.1 we 
discuss separation of Metric Inequalities, in subsection 3.2 we address the prob- 
lem of finding violated Metric Inequalities with minimum support. In 3.4 and 
3.5 we describe, respectively, separation of Metric Inequalities with p £ {0, 
and separation of Mod-/c Cuts. 

In section 3.6 we describe the shrinking procedure that, for a given metric p, 
allows computing p(p, G, d) on a reduced size graph. 
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3.1 Separation of Metric Inequalities 

It is well known that the separation problem for Metric Inequalities amounts 
to checking whether there exists a feasible multicommodity flow with respect 
to a given capacity vector x G R)_ . The feasibility of a multicommodity flow 
problem can be tested by solving the following Linear Programming (auxiliary) 
problem Sepl(G, d ): 



max a 



E /« 


- £ //« = —d°ot, ieV,oeO 


(6) 


j-.ijGE 


j.jieE 




E </« + 


fji)<Xij, ij£E 


(7) 


oGO 






fij > o, 


oeO,ij G E 





If the above problem admits an optimal solution of value a > 1 then the mul- 
ticommodity flow is feasible and no Metric Inequality is violated. If, conversely, 
the optimal solution has a value a < 1 then the multicommodity flow problem 
is infeasible and the optimal dual solution of the auxiliary problem provides a 
violated Metric Inequality. The above auxiliary problem is the well known Max- 
imum Concurrent Flow Problem (see Giinluk [26], Bienstock and Raskina [13] 
and Fleischer [24]). Its dual has the following form: 

min E %ij l^ij 

ijeE 

( fij ) C — + dij > o G O , ij G E 

ifji) £j ~ P-i + fj’ij ^0, O G O , ij £ E 

(a) J2Y, d ° e ° = 1 ( 8 ) 

oeo iev 

Pij > 0, ij G E 

Denoted by (p,,£) its optimal solution and assuming that a < 1, we have that 
xji < 1 and, consequently, that xfi < £d. This implies that fix > id is the sought 
for a maximally violated Metric Inequality. 

The following result proves that fi G Met(G). 

Lemma 1. Let (fi,£) be the optimal solution of the separation oracle Sepl{G, d). 
Then fi G Met(G). 

Proof. For any hk G E, let p{Pffff n ) denote the length of the shortest path 
between h and k according to the weights p. Suppose that fi £ Met(G). Then 
there exist ij G E such that fi t j > fi(Pij). But since 1° is the length of the 
shortest path from o to i, we can decrease fi t3 while keeping triangle inequalities 
satisfied. We get a better feasible solution and consequently (fi, £) is not optimal, 
contradiction. 
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3.2 Separation of Strong Metric Inequalities 

Estimating the quality of cutting planes is a major issue in Computational In- 
teger Programming as, in general, violation cannot be considered as a good cut 
selection criterion. Bienstock and Bley [9] observe that looking for minimal sup- 
port violated Cut Inequalities may lead to more effective cutting planes and 
suggest to perturb the fractional solution by a tiny value to get minimal support 
violated cutting planes. 

In what follows we generalize such considerations by measuring the depth of 
a Metric Inequality as the ratio between violation and the size of the support 

of the inequality, defined as E ^ . We refer to Metric Inequalities minimizing 

ijeE 

the ratio: 

E -EE 

ijeE oeQiev ^ 

E Pij 

ijeE 

as Strong Metric Inequalities. The separation oracle for Strong Metric Inequali- 
ties is: 



E XijfMj -EE 

. ijeE oeOiev 

mm — 

E tnj 

ijeE 

e° - i° + iMj >0, o G O, ij G E (10) 

I°j — T jiij > 0, o G O, ij G E 
Pij > 0, ij G E 

We prove that the separation oracle for Strong Metric Inequalities reduces 
to the Capacity Reduction Problem Sep2(G, d): 

max (3 

E fii- E = oeO,ieV 

j-ijeE j.jieE 

E (f°J + //*) ^ ~ P' C e E 

oeo 

f°> 0, oeO,ijeE 



If the above problem admits an optimal solution of value /3 > 0 then the mul- 
ticommodity flow is feasible and no Metric Inequality is violated. If, conversely, 
the optimal solution has a value (3 < 0 then, as in the previous case, the optimal 
dual solution will provide us with a violated constraint. The dual of Sep2(G, d), 
say Sep2dual(G,d) is: 
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min E X%j f-^ij 

ijeE oeOieV 

£° - e° + Hij >0, o G O, ij G E (11) 

H- (iij > 0, o G O, ij G E 

E = 1 

ij€E 

Hij > 0, ij G E 

Denoted by (/t, Z) its optimal solution and assuming that f3 < 0, we have that 
xfi — d£ < 0. This implies that fix > £d is the sought for maximally violated 
Metric Inequality. 



Theorem 3. The separation oracle (11), provides a Metric Inequality minimiz- 
ing the ratio (9). 



Proof. Consider the fractional linear program (10). By letting 1/ Hij 

ijeE 



we can write the equivalent problem: 



= t, 



min E tXij flij EE td°£° 

ijGE o&OieV 

t£° — t£° + tHij > 0, o G 0,ij G E 

t£° — t£° + tHij >0, o G 0,ij G E 

t E IJ P = 1 

ij£E 

Hij > 0, ij G E 

Then, by letting tHij = Wp and t£° = v°, we get: 



min XijWij - Y,Y< v ° d ° 
ijeE o&OieV 

v° — v° + Wij >0, oGO, ij G E 

v< j T >0, oGO, ij G E 

E Wi i = 1 

ij&E 

Wij > 0, ij G E 
t > 0 



So we have obtained Sep2dual(G , d) by linearizing the fractional linear pro- 
gramming problem (10). 

In table 1 we show for a set of test instances that separation of Strong Metric 
Inequalities yields a strong reduction in the number of violated cutting planes 
used to provide a lower bound. 
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Name 


LB 


jf 1 Metric 


#Strong 


Sun.trl 


2753.39 


402 


135 


Sun.tr2 


2791.53 


480 


175 


Sun.tr3 


3067.96 


473 


191 


Sun.tr4 


2829.51 


456 


163 


Sun.tr5 


2391.3 


412 


188 


Sun.tr6 


2997.64 


450 


175 



Table 1 . Comparison of separation oracles Sepl(G,d) and Sep2(G,d). 



3.3 A Lifting Procedure for Metric Inequalities 

Working on the fractional support S = {ij G E : Xij > 0} is crucial to make 
separation algorithms more efficient. 

Let Gs(V, S) be the fractional support graph. Here we present a lifting pro- 
cedure that extends a Metric Inequality defined on Gs to a globally valid Metric 
Inequality defined on G. 

Let {ps,i) be the optimal dual variables of the separation problem. The 
vector jig can be extended to the metric ji defined on G by letting: 

Hij = max{|£° — i °\, o G O }, ij G E\S 

The lifting procedure is exact for Sepl(G, d), i.e. it guarantees that the Metric 
Inequality ps x S > i d found on Gs and the globally valid Metric Inequality 
HS X S > id, obtained by the lifting procedure, have the same amount of violation. 
For the Strong Metric Inequalities the lifting procedure ensures that a violated 
Metric Inequality exists on G if and only if a violated Metric Inequality exists on 
Gs, but there is no guarantee that fix > id has the same ratio (9) as /. isxs > id. 

3.4 Separation of {0, l}-Rounded Metric Inequalities 

For some classes of instances it turned out to be effective formulating the sepa- 
ration problem of {0, 1 {-Rounded Metric Inequalities as a Mixed-Integer Linear 
Programming problem: 

min Xij Hij — z 

ij&E 

i°-ij+ jlij 0, O G 
ij ~ i° + i l ij >0, O G 

^ < EE d “ l ”+ 1 

oeo iev 
Hij G {0,1}, 

Z G ^-|- 

Constraint (12) enforces z = \id ]. If the optimal objective value is negative, we 
have found a violated {0, l}-Rounded Metric Inequality n x > z - 



O, ij G E 
O, ij G E 



ij G E 



(12) 
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We observe that {0, l}-Rounded Metric Inequalities include Cut Inequalities, 
which are Tight Metric Inequalities, as observed in Section 2. 

3.5 Mod-fc Cuts 

Given the integer polyhedron Pj = conv{ x £ Z n : Ax > 6}, where A £ Z mxn 
and b £ Z m , a Chvatal-Gomory cut is a valid inequality for Pj of the type 
A Ax > [A6] for some A £ such that XA £ Z". 

A Chvatal-Gomory Cut is a mod-k cut if A* £ {0, 1/fc, . . . , (k — 1 )/k} for 
some integer k > 2 ([14], [15], [29]). Every non dominated Chvatal-Gomory Cut 
is a Mod-fc Cut for some k. For any given x* £ P, the vector s* = Ax* — b is 
called slack vector. A cut is said to be totally tight at x* if s*Xj = 0 holds for 
j = 1, • • • ) m. 

The separation of a violated Chvatal-Gomory Cut was shown to be strongly 
NP-hard [23] . The same happens for Mod-fc cuts [14] , but to find a totally tight 
Mod-fc Cut whose violation is maximal, or prove that none exists, can be done in 
polynomial time, provided that the prime factorization of fc is known, by solving 
a congruence system in GP(fc), as described in [15]. 

If inequalities in the matrix A are Metric Inequalities, that is they have left- 
hand-side coefficients that define a metric, any inequality obtained by combining 
them is still a Metric Inequality, with a possibly larger r.lr.s. It follows that 
Chvatal-Gomory Cuts and Mod-fc Cuts can be seen as an intermediate step 
towards Tight Metric Inequalities. They also include, as special case, some Tight 
Metric Inequalities, like the Three-Partition Inequalities [32] . 

3.6 Computing p(p, G, d): a Shrinking Procedure 

Here we address step ii ) of the separation heuristic for Tight Metric Inequalities. 
Any valid inequality px > b , generated as described in sections 3. 1-3. 5, can be 
tightened by replacing the r.lr.s. b with p(p,G,d). 

To compute p(p, G , d) efficiently, we introduce a shrinking procedure 
that builds a reduced size problem NL(ph,Gh,dh) with the property that 
p(ph, G hl d h ) = p(p, G, d). 

Let ij £ E such that p t] = 0 and let G h (V h ,E h ) be the graph obtained by 
shrinking the nodes i and j into the supernode h, thus V h = V\{z, j}U{fc}. The 
edge set E h is defined by merging all the pair of edges vi, vj £ E into vh £ E h , 
for each v £ V \ {■ i,j }. 

Let p h denote the mapping of p onto G h . We observe that, due to the triangle 
inequalities p v j < p V i + p v j, p V i < p v j + p u j and pij = 0, we get p„ = p V j. For 
any v £ V h \ {h} we set p^ h = p V i = p v j and p^ v = p uv , for each edge of the 
set {uv £ E : u,v £ V \ {i, j}}. The following result is due to Lomonosov [30]: 

Theorem 4. Let p be a metric. The vector p h is a metric if Pij = 0. 

Accordingly, let D h be graph obtained by shrinking the nodes i and j in the 
demand graph D and let dh be the mapping of d onto D h . We have df JV = d uv 
for any u, v £ V h \{h} 1 d% h = d vi + d vj and dfc v = d iv + d jv . 
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Fig. 1. The shrinking procedure. 



Let p{p h ,G h ,d h ) denote the optimal solution of the reduced problem 
NL(p h ,G h ,d h ). The following property holds 

Theorem 5. Let p be a metric on G. Let ij £ E such that pij = 0 and let 
p h ,g h ,d h be defined as above. Then p(p,G,d) = p(p h ,G h ,d h ). 

Proof. Suppose that p(p h ,G h ,d h ) < p(p,G,d) and let x h be the optimal so- 
lution corresponding to p(p h ,G h 7 d h ). We can construct a feasible solution x 
for the original problem NL(p,G,d) as follows. Let Xij = [ 'fT d uv ~\, let 

uv(zH 

x uv — Xuv> uv € E and let x V i + x V j = x!f h . We know that px = p h x h = 
p(p h ,G h ,d h ) < p(p,G,d), contradiction because p{p,G,d) is the minimum. 

Now suppose that p(p h ,G h , d h ) > p{g, G, d) and let x be the optimal solution 
of p(p, G, d). Let x h be the mapping of x onto G h . Then p h x h = px < p{p, G, d), 
contradiction because p(p h ,G h ,d h ) is the minimum. 

Applying the shrinking procedure iteratively we build a reduced size graph. 
If the size of the reduced graph is less than a given threshold (< 8 nodes in our 
experience), we can efficiently compute p{p, G, d). 



4 The Cutting Plane Algorithm 

A cutting plane algorithm has been implemented to validate the effectiveness of 
the proposed separation procedures. 

The skeleton of the algorithm is Cplex 8.1 that provides search tree and cut 
pool management functionalities. We adopt the default branching strategy of 
Cplex 8.1. Variable selection is performed by strong branching. 



4.1 Preprocessing 

A simple operation is used to delete redundant variables: let c(Pij ) be the length 
of the shortest path between i and j according to the objective function costs. 
We delete the edge ij if > c(Py). 
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Let us consider the optimal solution of the LP-relaxation, whose value will 
be denoted by LB. For any nonbasic variable we can fix bounds by reduced cost 
fixing. Let LB be the value of the LP-relaxation and let UB be the value of 
the best feasible solution found so far. Let r.-,j be the reduced cost of variable 
x,j . For each nonbasic variable x.jj we compute S = and we fix u ij = 

min(uij,lij + 6) and lij = max(lij,Uij — (5). 

4.2 Initial LP-relaxation 

The initial LP-relaxation includes the Degree Inequalities x(6(i)) > |"d(<5(i))] , 
i G V. 

4.3 Separation Strategy 

At the root node of the Branclr-and-Cut tree we run separation routines in the 
following priority order: Rounded Strong Metric Inequalities, Mod-fc cuts. Once 
a violated cut has been found, we apply the shrinking procedure of section 3.6 
to get a Tight Metric Inequality. 

At the deeper nodes we look only for violated Rounded Strong Metric In- 
equalities to guarantee feasibility. 

Given the Metric Inequality fix > £d, we transform it into the Rounded 
Metric Inequality \\gf\x > \X£d] by letting A = min : ij G E}, as suggested 
in [10]. 

4.4 Computing Upper Bounds 

Let X = {x i: j : ij G E} be the optimal (fractional) solution of the LP-relaxation 
of the Capacity Formulation, i.e. X satisfies all the Metric Inequalities. A trivial 
upper bound can be computed by rounding-up the Xij, for ij G E. 

The trivial upper bound can be improved by a LP-based heuristic. We impose 
the bounds [xij\ and [Sjj], for each ij G E and we solve a restricted problem 
with the new bounds, yielding a minimum cost round-up of the Xij that preserves 
feasibility. 

Alternatively, the trival upper bound can be improved by a re-routing heuris- 
tic that attempts to round-down some of the fractional Xij. 

For each hk G E, let Shk = \xhk] ~~ Xhk be the slack of hk. 

For any candidate edge ij G E, the heuristic attempts to re-route the amount 
of flow fij = Xij — [Xij\ along a path Lfj from i to j formed by edges with 
Shk > fij,hk G Pij. If such a path exists, we set Xij = [xij\, yielding a new 
feasible fractional solution that provides a better upper bound. We found useful 
to sort the candidate edges ij in non-increasing order according to the value 
Cij ( \Xjj ] — Xij ) . 

At every node of the Branch-and-Cut tree we run a plunging heuristic. Let 
A be a fractional solution satisfying Metric Inequalities. If X is near-integer, i.e. 
if Xij < \xij\ + £ V Xij > \xij] — £, for ij G E, then we plung a 'ij to the the 
closest bound to get an integer solution Y and if Y is feasible we can update the 
current upper bound. 
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5 Computational Results 

In this section we report on the computational experience with the cutting plane 
algorithm. We ran the experiments on a Pentium IV 1.7 Ghz Personal Computer, 
with 512 Mb RAM, using Cplex 8.1 as LP solver and as Branch-and-Cut frame- 
work. 

The test bed is organized into two sets. The first set includes the well-known 
Norwegian instances, whose asymmetric version is also tested. The second set 
includes instances derived as backbone networks from two real-life instances 
introduced by Baralrona [4]. 

The format of the tables is the following. Name is the name of the in- 
stance, LB plow is the optimal value of the LP-relaxation of the Flow Formu- 
lation, LBcuts is the lower bound after the addition of cutting planes, BLB 
is the best lower bound yielded by the Branch-and-Cut, UBneu is the upper 
bound produced either by the LP-based or by the re-routing heuristic, BUB 
is best upper bound produced by the Branch-and-Cut (we use boldface to de- 
note the upper bounds proved to be optimal), Gap is the final gap, computed 
as 100 * ( BUB — BLB) / BU B), #RMet and #Mod-/c are respectively the total 
number of Rounded Strong Metric Inequalities and the total number of Mod-/c 
cuts added to the formulation. Finally, columns Nodes and Time report respec- 
tively on the number of nodes of the Branch-and-Cut tree and on total CPU 
time (secs.) spent in computation. 



5.1 Norwegian Instances 

Sun networks, with 27 nodes and 51 edges, were introduced in [21] and are 
known in the literature as the Norwegian networks. For each network we have 
tested the two different traffic patterns proposed in [10]: Sun.tr , 67 demand 
pairs with magnitude in the interval [0.1, 1.2] and Sun. dense, 702 demand pairs 
with magnitude in the interval [1, 2]. For the instances in this class the LP-based 
upper bound heuristic, whose results are reported in the column UBneu, turned 
out to be very effective, both in the symmetric and in the asymmetric case. 



Name 


LB Flow 


LBcuts 


BLB 


U BHeu 


BUB 


Gap% 


#RM 


#Mod-fc 


Nodes 


Time 


SUN.trl 


2738.78 


2822.52 


2825.44 


2830.74 


2825.44 


0.00 


138 


34 


31 


17 


SUN.tr2 


2773.72 


2862.19 


2869.92 


3150.98 


2865.30 


0.00 


171 


31 


19 


78 


SU N.tr3 


3038.38 


3139.19 


3153.13 


3153.67 


3153.13 


0.00 


160 


42 


87 


401 


SUN.trl 


2824.11 


2891.19 


2894.41 


2895.9 


2894.41 


0.00 


121 


32 


22 


632 


SU N.tr5 


2375.11 


2460.47 


2472.10 


2473.34 


2472.10 


0.00 


171 


142 


91 


58 


SU N.tr6 


2979.25 


3075.32 


3086.67 


3090.34 


3086.67 


0.00 


206 


23 


562 


168 



Table 2. Computational results for the SUN.tr instances 
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Name 


LB Flow 


LBcuts 


BLB 


U B Heu 


BUB 


Gap% 


#RM 


#Mod -k 


Nodes 




SUN.densel 


29594 


29650.9 


29651.3 


29651.3 


29651.3 




188 


33 


5 


mm\ 


SU N.dense2 


29583.9 


29646.6 


29653.6 


29654.2 


29653.6 


0.00 


282 


30 


663 


Eltfl 


SUN.dense3 


98575.5 


98639.8 


3153.1 


98640.5 


98640.5 




186 


30 


10 




SUN.denseA 


98347.9 


98406.0 


2894.4 


98414.4 


98410 


0.00 


209 


41 


130 




SU N.denseS 


59115 


59163.2 


59166.8 


59167.4 






217 


16 






SUN.denseG 


58923.4 


58985.6 


58987.8 


58990.2 






183 


30 




295 | 



Table 3. Computational results for the SUN. dense instances 



Asymmetric Norwegian Instances Asymmetric instances, denoted by 
Sun.tr and Sun.dense, are derived, respectively, from Sun.tr and Sun. dense 
by replacing every edge ij £ E with two directed arcs ij,ji with costs = Cji . 
Previous results for these instances are reported in [10] and, as far as our knowl- 
edge, none of these instances, with the only exception of Sun.trA , have been 
solved to exact optimality before. 

For the Sun.tr instances, Mod-fc cuts turned out to be ineffective and were 
replaced with the {0, l}-Rounded Metric Inequalities, whose separation is per- 
formed using the MIP problem of section 3.4. Let #01 RM denote the number 
of {0, l}-Rounded Metric Inequalities added in the search tree. 

In Table 4 we report computational results for the Sun.tr instances, that 
were solved to optimality. 



Name 


LBpi ow 


LBcuts 


BLB 


UBHeu 


BUB 


Gap% 


#RM 


#Mod-fc 


Nodes 


Time 


Sun.trl 


2738.78 


2962.42 


2962.42 


2979.49 


2962.42 


0.00 


3029 


161 


5914 


4523 


Sun.tr2 


2773.72 


2976.24 


2976.24 


2980.71 


2976.24 


0.00 


1684 


94 


1928 


1733 


Sun.tr 3 


3038.38 


3242.78 


3242.78 


3257.22 


3242.78 


0.00 


1406 


73 


426 


1880 


Sun.tr 4 


2824.11 


2978.90 


2978.90 


2978.9 


2978.90 


0.00 


506 


132 


46 


7103 


Sun.tr5 


2375.11 


2595.00 


2595.00 


2609.87 


2595.00 


0.00 


2033 


113 


1591 


8947 


Sun.trQ 


2979.25 


3196.96 


3196.96 


3196.96 


3196.96 


0.00 


8539 


108 


4487 


20464 



Table 4. Computational results for the Sun.tr qiustances 



Table 5 reports on results for Sun.dense instances. Due to the larger number 
of demand pairs and of source nodes, separation of Strong Metric Inequalities 
becomes more difficult. For such instances we were not able to prove exact op- 
timality, but we yielded small gaps and improved the results reported in [10]. 

5.2 US and German Instances 

In [4] Baralrona introduced two real-life instances defined on complete graphs, 
that were kindly made available to us by the author: BarAl, a US instance with 
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Name 


LBfIow 


LBcuts 


BLB 


UBneu 


BUB 


Gap% 


#RM 


#Mod-fc 


Nodes 


Time 


Sun.dense 1 


29594.0 


29721.8 


29732.2 


30641.0 


30265.1 


0.29 


8188 


23 


4487 


18000 


Sun.dense2 


29583.9 


29715.8 


29724.1 


30647.1 


30219.9 


0.22 


13626 


46 


1773 


18000 


Sun.dense 3 


98575.5 


98697.0 






99329.7 


0.06 


13588 


25 


1850 


18000 


Sun.dense 4 


98347.9 


98466.5 






99092.4 


0.56 


7789 


- 


555 


18000 


Sun.dense 5 


59115.5 


59228.8 


59244.2 


60157.9 


59847.5 


0.19 


14649 


26 


1687 


18000 


Sun.dense 6 


58923.4 


59038.4 


59051.2 


60014.1 


59667.5 


0.08 


11871 


25 


2627 


18000 



Table 5. Computational results for the Sun.dense instances. 



41 nodes and 646 demand pairs and Sar64, a German instance with 64 nodes 
and 2016 demand pairs. ^From each of these two instances a subset of nodes 
was selected to form two backbone networks with, respectively, 13 and 9 nodes. 
We generated more backbone networks, defined by the subgraph induced by 
the nodes with largest demands, that will be made public available through a 
website. 

Nevertheless we remark that BarAl and _Bar64 are still out of the scope of 
our algorithm as we were not able to solve the LP-relaxation. In table 1 we 
report on computational results for the backbone instances. We set a time limit 
of 10800 seconds for the search tree. Here the LP-based heuristic turned out to 
be ineffective and we used the re-routing heuristic, whose results are reported in 
the column UBneu- 



Name 


LB Flow 


LBcuts 


BLB 


U BHeu 


BUB 


Gap% 


#RM 


#Mod-fc 


Nodes 


Time 


BarAl/9 


155205 


174517 


176073 


176385 


176073 




90 


9 


96 


37 


Baril/12 


251702 


286060 


291188 


294196 


291188 


0.00 


524 


35 


96 


480 


BarAA/19 


307044 


352559 


359237 


363945 


363945 


2.2 


KEfifil 


74 


685 


12859 


BarAA/2\ 


316817 


365739 


368338 


474059 


474059 


22.3 


4297 


121 


2236 


219 


OO 

c* 

s- 

e 

cq_ 


383692 


438536 


438536 


594692 


594692 


26.2 


21238 


194 


2236 


118500 


Bar64/13 


686817 


889686 


961550 


977750 


961550 




2184 


27 


4711 


3110 


Bar64/14 


705518 


937792 




1011980 






696 


38 


1443 


1008 


Bar64/17 


764169 


995694 


[jjj] 


1099800 






2332 


52 


5914 


372 


Bar64/19 


807992 


1070610 


1134030 


1623100 


1623100 


30.1 


2886 


71 


2882 




Bar64/26 


951923 


1262090 


1301920 


2089430 


2089430 


37.7 


2975 


99 


718 


mm 


Bar64/32 


1076730 


1477130 


1495760 


2406900 


2406900 


37.9 


5693 


125 


137 


12871 | 



Table 6. Computational results for the BarQA/ XX and BarAl/XX instances 
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Abstract. In this paper we use facets of the convex hull of mixed- integer 
sets with two and three variables to derive valid inequalities for integer 
sets defined by a single equation. These inequalities also define facets 
of the master cyclic group polyhedron of Gomory. In particular, our 
inequalities generalize the 2slope facets of Araoz, Gomory, Johnson and 
Evans (2003). In addition, they dominate the strong fractional cuts of 
Letchford and Lodi (2002). 



1 Introduction 



An important technique in generating cutting planes for integer programs is to 
work with a single implied constraint, and to derive valid inequalities for the 
non-negative integral vectors satisfying the constraint. Suppose we express the 
constraint in the form J a j x j + U — b, and let 

Y + = ja; G Z^ J \ y G Z : a-jXj + y = b, x > 0, y > 0 j. 



Valid inequalities for Y + can be used as cutting planes for the original integer 
program. One way to generate valid inequalities for Y + is to use facets of the 
master cyclic group polyhedron: 



P(n, r) = convex 



G Z 



n— 1 



n— 1 

£■ 

i= 1 



;= r(modn), w > oj 



(1) 



where n,r G Z and n > r > 0 and a= 6(modn) means a — b is divisible by 
n. For a set S C R n , conv(S) denotes the convex hull of vectors in S. P(n,r) 
is closely related to a relaxation of Y + obtained by relaxing the non-negativity 
requirement on the y variable: 

Y = | a: G Z^ J \ y G Z : aj Xj + y = b , x > oj. (2) 

Facets of P(n,r ) yield inequalities valid for Y , which are also valid for Y + . 
The polyhedral structure of P(n,r) was initially studied in Gomory [11] and in 



D. Bienstock and G. Nemhauser (Eds.): IPCO 2004, LNCS 3064, pp. 33—45, 2004. 
(c) Springer- Verlag Berlin Heidelberg 2004 




34 



S. Dash and O. Giinliik 



Gomory and Johnson [12,13]. More recent references include [15], [3], and [14]. 
Also see [1] for an introduction. 

In this paper we study facets of P(n, r) and valid inequalities for Y that can 
be derived from the following simple polyhedral sets: 

Q 1 = {x £ R, y £ Z : x + y > b , x > 0} 

and 

Q 2 = {x £ R, y,z £ Z : x + ay + z > (3, x, y > 0 } , 

and their relaxations obtained by allowing the integral variables in these models 
to be (l/f)-integral for t £ Z. Our work is partly motivated by the work of 
Marchand and Wolsey[18], who derive the well-known Gomory mixed-integer cut 
(GMIC) using Q 1 . GMIC can also be derived from a facet of the master cyclic 
group polyhedron [11], and our interest in P(n, r) is partly due to the practical 
importance of the GMIC, which is now routinely used in MIP software. 

As discussed in Gomory and Johnson [12,13], there is a one-to-one corre- 
spondence between the facets of P(n, r) and the facets of the mixed-integer set 
obtained by including continuous variables in P(n,r). From this, one can show 
that facets of P(n,r) yield cutting planes for general mixed-integer programs. 

Most of the proofs have been omitted in this paper due to space restrictions, 
but they can be found in [8]. 



2 Scaled MIR Inequalities 



The simple mixed-integer rounding (MIR) procedure provides a unifying frame- 
work to derive valid inequalities for mixed-integer sets. See [18] for examples. 
Given a valid inequality x + y > b for a mixed-integer set X C { (x, y) : x £ 
P+, y £ Z}, the MIR inequality x > r(\b ] — y) where r = b — [b\ is also valid 
for X, see [19]. 

The MIR inequality is facet defining for the set conv(Q 1 ), and 

conv^Q 1 ) = {x, y £ R : x + y > 6, x > r([6] — y), x > 0}. 

We emphasize that the variable y in Q 1 can be negative. For any t £ Z we next 
define a relaxation of Y by letting the y variable take on 1/f-integral values: 

Y l = {x £ Z^ J \ty £ Z : aj Xj + y = b , x > 0}. 



Multiplying the equation by t and substituting z 1 = ty , we observe that the 
resulting set has the the same form as Y . Define c* to be tc — [tc\ for c £ R and 
t £ Z . If tb is not integral, we can use the MIR inequality, as in the derivation 
of the GMIC in Marchand and Wolsey[18], to show that 



^ . P 

j&J : 



E 



164 



o‘.>6‘ 

j — 



i z3. 

1 — 6 * 



b > 



(3) 
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is valid for Y* . Therefore, inequality (3), which we call the t-scaled MIR inequality 
or the t-MIR cut, is valid for Y . 

When applied to a row of the simplex tableau, the 1-scaled MIR inequal- 
ity gives the well-known Gomory mixed-integer cut. When t > 1, the resulting 
inequality is referred to as the k-cut in Cornuejols, Li and Vandenbussche [7]. 
We next define the t-scaled MIR function and show that there are only a small 
number of t-scaled MIR inequalities for Y . 

Definition 1 . For a non-zero integer t, the t-scaled MIR function with param- 
eter b is defined as: 



f t,b (v] = f v*/ b* if it* Kb*, 

3 \ (1 — v*)/(l — b 1 ) if v* > b*, 

where v* = tv — [tuj and b* = tb — [tb\ . 

Using this definition, the t-scaled MIR inequality (3) becomes 

J2f t,b ( a o) x o ^ 1 - 

j£J 



Lemma 1. Assume the inequality defining Y is rational, and let n be the small- 
est integer such that nb £ Z and naj £ Z for all j £ J . Then, there are at most 
|_n/2j distinct t-MIR cuts for Y. In particular, if n > t > \n/2\ then the t-MIR 
cut is the same as the (n — t)-MIR cut. □ 

We note that Example 1 in Cornuejols, Li and Vandenbussche [7] deals with the 
set Y = { x £ Z 2 , y £ Z : 1.4aq + 0.1:r2 + y = 4.3, x > 0} which has at most 
5 distinct t-MIR cuts due to Lemma 1 because n = 10. This explains why the 
scaled cuts for f and 10 — t turn out to be identical in [7]. 

It is not difficult to see that the equation defining P(n,r) can be written in 
the same form as Y and therefore t-MIR cuts can be applied to P(n, r ). Gomory 
[11] showed that the 1-MIR cut defines a facet of P(n,r). Combining results on 
automorphisms and homomorphisms in [11], we can prove the following result. 

Theorem 1. For every integer t = 1, 2, ... , |_rz/2j , such that that tr is not an 
integral multiple of n, the t-MIR. cut defines a facet of P(n,r). □ 

Gomory, Johnson and Evans [15] present results based on a shooting ex- 
periment of Gomory and identify the “important” facets of the cyclic group 
polylreclra for small n (i.e., n < 20). In particular, they present coefficients of 
13 important facet defining inequalities for P(10,9), P(20, 19), P(15,14), and 
P(13, 12). Intriguingly, all of these facets are scaled MIR facets. In Table 2, we 
list these inequalities with the corresponding master polyhedra, and their scaling 
parameter. 
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Relative importance of facet 


Polyhedron 


Reference in [15] 


1 


2 


3 


4 


P(10,9) 


Table 2 


5-MIR 


2- MIR 


4-MIR 


- 


P(20,19) 


Table 3 


10- MIR 


5-MIR 


4-MIR 


- 


P(15,14) 


Figures 5-8 


5-MIR 


3-MIR 


6- MIR 


1-MIR 


P(13,12) 


Figures 9-11 


1-MIR 


2- MIR 


6-MIR 


- 



Table 1 . Important group facets are scaled MIR 



3 A Simple Polyhedral Set 

In this section we study simple three variable mixed-integer sets. We first look 
at the set 

Q 2+ = {x £ R, y, z £ Z : x + ay + z > j3, x, y, z > 0}, 

when a,/3 £ R and satisfy 1 > (3 > a > 0, and [/3/a] > (3 /a. Note that variable 
z is required to be non-negative. In a recent study Agra and Constantino [2] 
study Q 2+ when (3 is an arbitrary positive number and give a polynomial-time 
algorithm that enumerates all of its facets. Under our restrictions on a and (3 , 
it is possible to describe the convex hull of Q 2+ explicitly. 

Lemma 2. The following inequalities are facet defining for conv(Q 2+ ), and 
together with the non-negativity requirements, give a complete description of 
conv{Q 2+ ): 

x + ay + /3z> (3, (4) 

(l/(f3-a[/3/a\))x + y+ \/3/a~\ z > \ /3/a] . (5) 

Proof, (validity) Note that inequality (4) can be obtained by treating x + ay as a 
continuous variable and writing the MIR inequality based on (x + ay) + z > (3. 
To obtain inequality (5), we start with inequality (4), divide it by a and relax 
the resulting inequality as follows: 

x/a + y + \ f3/a) z > f3/a (6) 

Writing the MIR inequality where x/a is treated as a continuous variable and 
y + \/3/a \ z is treated as an integer variable gives inequality (5). □ 

Let Q 2 be the relaxation of Q 2+ obtained by allowing the 2 variable to assume 
negative values: 

Q 2 = {x £ R, y , z £ Z : x + ay + z > (3, x,y > 0} 

and remember that (3,a £ R satisfy 1 > (3 > a > 0, and | (3 /a] > (3 /a. Even 
though (3 is required to be less than 1 in Q 2 , the fact that 2 can take on negative 
values makes the set fairly general. We next show that (4) and (5) are facet 
defining for Q 2 under mild conditions. They are not necessarily sufficient to 
describe its convex hull. 
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Lemma 3. Inequality (4) is facet defining for conv(Q 2 ). In addition, inequal- 
ity (5) is facet defining for Q 2 ifl/a> \(d/a\. 

Proof, (validity) Inequality (4) is valid as it is derived in the proof of Lemma 2 
without assuming that 2 is non-negative. To see that inequality (5) is valid, 
notice that the following inequalities are valid for Q 2 : 

(l/a)x + y + (/ 3/a)z > (d/a, 

(l/a)x + y + (1 /a)z > (d/a, 

and therefore, for any 7 £ R satisfying 1/a > 7 > (d/a, the following inequality 

(1 /a)x + y + 'yz > (d/a 

is also valid as it can be obtained as a convex combination of valid inequalities. 
As 1/a > \(d/a~\ by assumption, inequality (6) is valid for Q 2 , and therefore so 
is inequality (5). □ 

When 1/a = [7 3/a], inequality (6) is the same as (l/a)x + y+(l/a)z > (d/a, 
and hence inequality (5) is an MIR inequality obtained after scaling x+ay+z > (d 
by 1/a. 

Inequality (5) leads to strong inequalities even when all variables are non- 
negative. Let 

Q 3 = {x £ R, y,z £ Z : x + ay + z > b, x, y,z > 0}, 

where b = [6J + b > 1 with b ^ Z and b > a > 0. Furthermore, assume that 
1/a > \b/a~\ > b/a. If we treat x + ay as a continuous variable and 2 as an 
integral variable and apply the MIR procedure, we obtain x + ay + bz > b\b] 
which can be relaxed to 

(1 /a)x + y+ b/a z > b\b]/a. (7) 

The MIR procedure can be applied again by treating (l/a)a; as a continuous 
variable and ( y + \b/a 1 \ z) as an integral variable. Another possibility is to view 
the initial inequality as x + ay +(z—[b\) > b and obtain x + ay + b(z—\b\) > b 
as the first MIR inequality which then leads to 

(1 /a)x + y+ b/a (z — \b\) > b/a, (8) 

via convex combinations, for the second MIR step. Inequality (8) is strictly 
stronger than inequality (7) as b [6] /a = b/a+b[b\/a < b/a+ \b/a~\ [6J . Therefore 
the MIR procedure based on inequality (8) gives a stronger valid inequality. We 
illustrate this in Example 1 and in Figure 1. 

Example 1. Consider the set Q 3 with a = 0.4 and b = 1.7 and note that 1/a = 

2.5 > | b/a] = [0.7/0. 4] = 2. In this case, inequality (7) becomes 2.5x + y + 2z > 

3.5 whereas inequality (8) is 2.5x+y+2z > 3.75. The second step MIR inequality 

based on inequality (7) is 2.5a: + 0.5 y + z > 2 which is weaker than the second 
step MIR inequality based on inequality (8): 2.5a: + 0.75 y + 1.5 z > 3. We also 
note this last inequality is facet defining for this instance. □ 
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z 




Fig. 1 . Different relaxations of MIRs 



In Figure 1, we plot four inequalities for the example above, by showing the 
boundaries of their feasible regions: the original inequality; the MIR inequality 
based on 2 being integral; inequality (7), shown by a dotted line; and inequal- 
ity (8), shown by a dashed line. It is only when we apply the second MIR step 
that we distinguish between y and x\ therefore we can depict the above inequal- 
ities in the plane, the horizontal axis representing 2.5a; + y and the vertical axis 
standing for z values. 

4 Two-Step MIR Inequalities for T and P(n,r ) 

In this section we apply inequality (5) to obtain valid inequalities for Y and 
facets of P(n , r). Let w = l a j\ x j + y~ IAJ an d % = a j ~ l a j\ , b = b — \ b\. 
The equality defining Y can now be written as 

Y aj Xj+w = b. (9) 

jeJ 

Note that w is integral provided that x and y are integral. 

Let a £ Rbe such that b > a > 0, and 1/a > \b/a] > b/a. Let { Ji, J2, J 3 , J4} 
be a partition of the index set J such that (i) if j £ J2, then cij < a, and (ii) if 
j € J3, then hj > a. We can now relax (9) to obtain 

Y % x i + Y a Xj + Y2 (a + ( aj — a)) Xj + Xj + w > b. 

ieJi jeJ 2 jeJ3 jeJi 

Rearranging the terms leads to 

( Y & 3 x i + Y “ a ) x i) + “( Y x i) + (^Z X 3 +W )-^ 

jeJi 3 &J 3 je J2UJ3 jeJ 4 
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which resembles the set Q 2 as the first term is non-negative, the second term 
gives positive integral multiples of a, and the last term is integral. As we chose 
a to satisfy 1/a > [6/a] , we can write the following valid inequality for Y based 
on inequality (5): 



E % x o 

jeJi 



E (“i _ a ) X J ^ 

j&Js 

b/a 



(b — a b/a J f 



E 

J&J2UJ3 



b/a ( 



E x i + w 

j&Ji 



Letting p = 6 — a[6/aj and r = [6/a], we can show that the strongest inequality 
of this form is obtained by partitioning J into J[ = [j € J : % < p}, J| = 
{j G J : p < a,j < a}, .// = {j G J : a < dj < p(r — 1) + a}, J| = {j G J : 
p(r — 1) + a > a,j}. 

If r > 2, we strengthen the previous inequality by using a stronger relaxation 
of inequality (9) as the starting inequality. We partition the index set J into 
2 t subsets { Ji, J\, , , . . . , .7/ } as follows: For a given index j G J, let kj = 

min] [dj /a] , r} — 1. Index j is put in the set J® if a : j — kja < p; it is put in the 
set Jj. if aj — kja > p. We now relax (9) to obtain 



E *3 x 3 + E ( E kaxj + (ha + ( dj — ka )) Xj'j + Xj + w > 6, 

ieJ o fc= 1 jeJl jeJ ° +1 ieJt 

which can be rewritten as 



( E % + E E ( a,j — ka) XjJ + a (E £ kxjj + 

jeJ o k =i je j ° +1 fc=?be^uj ° +1 

( E x i + w ) - b - 



Applying inequality (5) and substituting for w leads to inequality ’Ybj^jljXj — 
pr y > pr [6], where 



7 j=pr [aj\ + 



( aj — ka) + k p if j G J° +1 , k = 0, . . . , r — 1 
fcp if j G J/, fc = l,...,r. 



We next define the two-step MIR function and formally state that the two- 
step MIR inequality is valid for Y . 



Definition 2. Let b,a G R be such that 6 > a > 0, and 1/a > [6/a] > 6/a. 
Define p = b — a\b/a\, and t = [6/a] . The two-step MIR function for a right- 
hand side b with parameter a is defined by 
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9 b ’ a (v) 



0(1 — pr) — k(v)(a — p) 
Pr{ 1 - b) 
k(v) + 1 — TV 
r(l - b) 



k(v)a < p 
k(v)a > p, 



where k(v) = min{ \v/oi\ , r} — 1, and v = v — |_uj and b = b — |_&| • 

Lemma 4. For any a £ R that satisfies the conditions in Definition 2, two-step 
MIR inequality for right-hand side b with parameter a 

J2g b ' a ( a j)xj >1 (10) 

j&J 

is valid for Y . □ 

Observe that g b,a (v) is a piecewise linear function with two distinct slopes, 
(1 — pr)/{pr{ 1 — b)) and —1/(1 — b). Also, 1 > g b,a (v) > 0, for all v £ R. In 
[13,14], Gomory and Johnson describe a family of piecewise linear “two-slope” 
functions, containing g b ’ a (v), which yields valid inequalities for Y. Lemma 4 
shows that, of the two-slope functions in [13], the ones of the form g b,a {v) can 
be derived from the the two-step MIR inequalities. 

It is possible to write t-scaled two-step MIR inequalities for Y by scaling 
the initial equality by t £ Z. For the sake of completeness, we next present the 
f-scaled two-step MIR inequalities for Y . Note that 1 > P > 0 even when t < 0. 

Lemma 5. Let t £ Z and a £ R be such that b* > a > 0, and 1/a > {b* /a ] > 
bf/a. Then, the t-scaled two-step MIR inequality j g tb,a (taj)xj > 1 is valid 
for Y. □ 

When applied to P(n, r), the two-step MIR inequalities derived in the previ- 
ous section yield a wide range of facets. In [3], the authors present a class of facets 
of P(n,r) which they call “2slope” facets. We will see that these 2slope facets 
are (n — l)-scaled (or (— l)-scaled) two-step MIR inequalities for appropriate 
choices of a. Initially we consider 1-scaled inequalities and present a result that 
generalizes Theorem 3.5 of [3[. We refer to the facets in Theorem 2 as general 
two-slope facets, in keeping with the notation in [3] . 

Theorem 2. Let A £ Z + be such that r > A > 0. The two-step MIR inequality 

n— 1 

Y,g r/n ’ A/n (i/n) Wi > 1 

2—1 

defines a facet of P(n,r) provided that n > A [r/Z\] > r. We call the corre- 
sponding facet the 1-scaled two-step MIR facet of Pin, r) with parameter A. 

We next show that t-scaled two-step MIR inequalities define facets of P(n, r) 
under some conditions. Let pf n (i) = ti mod n for integers i. 
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Theorem 3. Let t £ Z be such that tr is not divisible by n. Let A £ Z + be such 
that A = tk mod n for some integer k. The t-scaled two-step MIR inequality 

n 

Y,g tr/n ’ A/n (tiMwi>i 

i—1 

defines a facet of P(n,r) provided that /jf n (r) > A > 0 7 and n > A \// n (r)/A\ > 

MnW- 

Based on this result, we can show the following result. 

Corollary 1. The 2slope facets of P(n,r) in [3, Theorem 3-4] with parameter 
d are (n — 1 )-scaled two-step MIR facets of P(n,r) with parameter A = n — d 
such that t = 2, where t — \ \(n — r)/ A~\ . 

Figure 2(a) shows a two-step MIR facet for P(10, 7) obtained by letting 
A = 4. Figure 2(b) shows a 3-scaled two-step MIR facet for P(10, 9) obtained 
by setting A = 4. In Figure 2, the marked points denote the coefficients of 
the corresponding facet-defining inequalities and the solid lines represent the 
underlying two-step MIR functions. 




Fig. 2. Two different two-step MIR facets. 



5 Extended Two-Step MIR Inequalities 

In this section we derive new valid inequalities by considering the limiting be- 
havior of the two-step MIR function g b,a . For a given right-hand side b , the set 
of admissible values of a that give valid inequalities for Y is not clear from Def- 
inition 2. Let b £ R be fixed and let T denote the set of values of a that satisfy 
the conditions of Definition 2. Let T = where T d is the set of values of 

a such that \b/a] = d. In other words, T d = (b/d, 1/d] fl ( b/d , b/(d — 1)). 

Notice that T d ^ 0 for all d > 2. Also, a = 1/d is an admissible choice if and 
only if 1/d < b/(d — 1). Let d be the largest integer strictly less than 1/(1 — b), 




42 



S. Dash and O. Giinliik 



that is, d = 



1/(1 - 6 ) 



— 1. Then, one can show that a = 1/d is an admissible 



choice in Lemma 4 only when 2 < d < d. Also, 



( (b/d, 1/d] if 2 <d<d 

?'={ „ „ ( 11 ) 

{(b/d,b/(d- 1)) ifd>d. 

Let £ = {b/d : d > d,d £ Z} and notice that = (0, 1 /d] \ £. Next, we 

will derive extended two-step MIR inequalities for Y, by considering the limiting 
behavior of the function g b,a when a tends to a point p in the set £. We will 
separately consider the case when a converges to p from below and from above, 
as the limits are different. We use e — > 0 + to denote that e takes only positive 
values as it tends to 0. 



Definition 3. Let b £ R be such that b > 0. Let d be an integer such that 



1/(1 — b) — 1. The extended two-step MIR function for a right-hand side 



d > 



b with parameter d is defined by 



v/b 



h b ’ d (v)= \ /(v) + l-(d+l)t 



(d+l)(l-6) 
where l(v) = min{ I dv/b I , d}. 



if v is an integral multiple of b/d and v < b, 



otherwise. 



Lemma 6. For any integer d satisfying the conditions of Definition 3, and for 
any v £ R, 

lim gh’W^-Dv) = h b ' d (v). 

e — * 0 + 



To complete the discussion on the limiting behavior of g b,a , we note that for 
any integer d > d, the two-step MIR function converges to the (1-scaled) MIR 
function when a in g bi<x converges to a point in £ from above. 

Let {ai} C R n be a sequence converging to the vector a, and let {bi} C R 
be a sequence of numbers converging to b, such that afx < bi is valid for a 
polyhedron P, for all i > 0. Then a T x < b is also valid for P. This fact and 
Lemma 6 imply the following result. 

Lemma 7. Let t be an arbitrary integer, and let d be a positive integer that 
satisfies the conditions in Definition 3 with b replaced by tb. Then the following 
inequality, called the t-scaled extended two-step MIR inequality, is valid for Y : 

htb ’ d ( ta i) x o ^ 1 



(12) 
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5.1 The Strong Fractional Cut of Letchford and Lodi 

Let Y + be the restriction of the set Y obtained by requiring the y variable to 
be non-negative, i.e., Y + = Y D {y > 0}. In [16], Letchford and Lodi present a 
valid inequality for Y + which they call the strong fractional cut. Their inequality 
dominates the so-called Gomory fractional cut 

22 Xi > b. (13) 

*e J 

In this section we will show that their inequality is implied by the extended 
two-step MIR inequalities. For convenience, we first present their main result in 
our notation. 

Theorem 4. (Letchford and Lodi [16]) Suppose 6 > 0 and let k > 1 be the 
unique integer such that jrjrf < 6 < ■%. Partition J into classes Qo, ■ ■ ■ ,Qk as 
follows. Let Qo = {i £ J : a, < 6} and, for p = 1 ,k, let Q p = {i £ J : 
6 + (p— 1)(1 — b)/k < di < b + p( 1 — b)/k}. Then the following strong fractional 
cut is valid for Y: 



22 &iXi + 22 H (di - p/(k + 1 ))xi > b. 
i&Q o p=l ieQp 



(14) 



It can be shown that k = 



as 



H I ^ max {°’ 



1/6 1 — 1, and that inequality (14) can be written 
k(di - b) 



ieJ 



1 — 6 



} 



> b. 



(15) 



We now show that the (-l)-scaled extended 2-step MIR inequality with pa- 
rameter k dominates inequality (15). First consider a relaxation of h b,d defined 
by: 



h b ’ d (v) = 



l(v) + 1 - (d + l)z 

(d + i)(i-6) 



where l(v) = min{ dv/b ,c?}, Vv £ R. 



When v is an integer multiple of b/d and v < b, we know that h b d (v ) = v/b, 
and therefore h b,d (v) < h b,d (v). For all other values of v, h b ' d (v) and h b,d (v) 
are identical. Now, multiply the equation defining the set Y + by —1 and apply 
h~ b,k , where k = 1/6 — 1, to get the valid inequality h~ b ’ k (—aj)xj > 1. 

As (— aj) — [~aj\ = 1 — dj and (—6) — |_ — = 1 — 6, this is: 



E 

ie J 



i {k, 



k(l—ai) 

1-6 



} + 1 — (k + 1)(1 — di) 



(k + 1)6 



Xi> 1 



By simple algebraic transformations, the above inequality is the same as inequal- 
ity inequality (15). 
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Recall d in Definition 3. As we are working with the equation defining Y + 
scaled by —1, k in (15) is actually equal to d. This implies the following result: 



Proposition 1. The strong fractional cut (If) of Theorem f is dominated by the 
(-l)-scaled extended two-step MIR inequality (12). Furthermore, inequality (If) 



is valid for all k > 



m 



- 1 . 



6 Concluding Remarks 

In this paper we derived the two-step MIR inequalities (10) and showed how 
these inequalities generalize both the 2slope facets of P(n, r) in [3] and the strong 
fractional cuts in [16]. In deriving these inequalities, we considered only one of 
the many facets of the set Q 2 , and applied this facet to a specific relaxation of 
Y. It is somewhat surprising that such a seemingly narrow procedure results in 
such general inequalities. 

It is relatively straightforward to generalize the two-step MIR inequalities 
when the equation defining the set Y contains continuous non-negative vari- 
ables in addition to integral variables. An important question is whether these 
inequalities will be useful for solving general mixed-integer programs, the way 
GMICs are. Based on the observation that the two-step MIR inequalities have a 
number of features in common with the GMIC, we are optimistic. In particular, 
we note they are both based on facets of simple mixed-integer sets, and they 
define important facets of low-dimensional cyclic group polylredra. 
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Abstract. In this paper we characterize the “price of anarchy”, i.e., 
the inefficiency between user and system optimal solutions, when costs 
are non-separable, asymmetric and nonlinear, generalizing earlier work 
that has addressed “the price of anarchy” under separable costs. This 
generalization models traffic equilibria, competitive multi-period pricing 
and competitive supply chains. The bounds established in this paper are 
tight and explicitly account for the degree of asymmetry and nonlinearity 
of the cost function. We introduce an alternate proof method for pro- 
viding bounds that uses ideas from semidefinite optimization. Finally, in 
the context of multi-period pricing our analysis establishes that user and 
system optimal solutions coincide. 

Keywords: System and User-Optimization, Traffic Equilibrium, Price 
of Anarchy. 



1 Introduction 

There has been an increasing literature in the recent years trying to quantify the 
inefficiency of Nash equilibrium problems (user-optimization) in non-cooperative 
games. The fact that there is not full efficiency in the system is well known both 
in the economics but also in the transportation literature (see [1], [8]). This in- 
efficiency of user-optimization was first quantified by Papadimitriou and Kout- 
soupias [15] in the context of a load balancing game. They coined the term “the 
price of anarchy” for characterizing the degree of efficiency loss. Subsequently, 
Roughgarden and Tardos [21] and Roughgarden [20] applied this idea to the clas- 
sical network equilibrium problem in transportation with arc cost functions that 
are separable in terms of the arc flows. They established worst case bounds for 
measuring this inefficiency for affine separable cost functions and subsequently 
for special classes of separable nonlinear ones (such as polynomials). It should 
be noted that Marcotte presented in [16], results on the “price of anarchy” for 
a bilevel network design model. Recently, Jolrari and Tsitsiklis [14] also studied 
this problem in the context of resource allocation between users sharing a com- 
mon resource. In their case the problem also reduces to one where each player 
has a separable payoff function. Correa, Schulz and Stier Moses [4] have also 
studied “the price of anarchy” in the context of transportation for capacitated 
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networks. The cost functions they consider are also separable. The paper by 
Chau and Sim [2] has recently considered the case of nonseparable, symmetric 
cost functions giving rise to the same bound as Rouglrgarden and Tardos [21]. 
Wardrop [24] was perhaps the first to state the equilibrium principles in the 
context of transportation. Dafermos and Sparrow [5] coined the terms “user- 
optimized” and “system-optimized” in order to distinguish between Nash equi- 
librium where users act unilaterally in their own self interest versus when users 
are forced to select the routes that optimize the total network efficiency. Smith 
[22] and Dafermos [7] recognized that this problem can be cast as a variational 
inequality. In [6] Dafermos considered how the decentralized “user-optimized” 
problem can become a centralized “system optimization” problem through the 
imposition of tolls. Recently, Hearn and co-authors [12], [13], have studied the 
problem of imposing tolls in order to induce a behavior to users so that their 
route choices are optimizing the overall system. They study a variety of criteria 
for imposing tolls. Finally, Cole, Dodis and Rouglrgarden [3] have also studied 
the notion of introducing tolls (taxes) in order to make the decentralized prob- 
lem efficient in a centralized manner. The review paper by Florian and Hearn 
[10], the book by Nagurney [17], and the references therein summarize the rele- 
vant literature in traffic equilibrium problems. Traffic equilibrium problems are 
typically modeled through variational inequalities. The books by Facchinei and 
Pang [9] summarize the developments in the area of variational inequalities. 
Nash equilibrium problems arise in a variety of settings and model competitive 
and non-cooperative behavior (see Section 2 for more details). In this paper we 
study the inefficiency of equilibrium by comparing how the presence of competi- 
tion affects the total profit in the system in a decentralized (user-optimized) ver- 
sus a centralized optimization (system-optimized) setting. We establish a bound 
on the ratio of the overall profit of the system in these two settings. This work is 
the first to consider non-separable, asymmetric cost functions and is important 
since, as we discuss in Section 2, it allows modeling more realistic situations. Fur- 
thermore, it allows a unifying framework which naturally extends results in the 
current literature. In particular, our contributions versus the existing literature 
are the following. 

1. We consider non-separable functions in the sense that cost functions also 
depend on the strategies of the competitors. Furthermore, cost functions 
can be asymmetric in the sense that different competitors’ strategies affect 
their cost functions differently. This generalization is important since the 
strategies of one’s competitors will influence his/her own cost in an asym- 
metric way. In particular, we introduce a measure of asymmetry (denoted 
by c 2 in Section 2) which quantifies the degree of asymmetry of the 
competitors’ cost functions. We establish that the ratio of the total cost in 
the system operating in a user-optimized setting versus the total cost in a 
system optimized setting is bounded by 

ifc!£2 
c 2 if c 2 > 2. 
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We illustrate how our results are a natural generalization of the bound es- 
tablished by [20], [21] since when the problem functions are separable and 
affine the bound becomes 4/3. Furthermore, the work by Chau and Sim [2] 
(where the functions are non-separable but symmetric) becomes a special 
case since the bound is still 4/3, since c 2 = 1. The results in the affine case 
allow the feasible region to be a non-convex set. 

2. When the cost function has no constant term then the bound becomes c 2 . 
An important implication of this result arises in the context of multi-period 
pricing. In this setting our analysis establishes that user and system optimal 
solutions coincide. 

3. We generalize our results to nonlinear functions. We introduce a measure 
which quantifies the degree of nonlinearity of the problem function (de- 
noted by A in Section 3). We establish that the bound naturally extends to 
involve the nonlinearity parameter A, i.e., 

J 4^3 i! ^<i 

[c 2 A 2 - 2{A-1) if c 2 > 

4. We establish that the bounds above are tight for affine as well as for some 
nonlinear problems. 

5. We introduce an alternative semidefinite optimization formulation for 
deriving these bounds. This approach does not require positive definiteness 
of the Jacobian matrix (i.e., it does not need to be invertible). Therefore, the 
solution does not need to be unique. We illustrate that this approach gives 
rise to the same bound when the Jacobian matrix is positive definite. 

2 A Bound for Affine and Asymmetric Cost Functions 

In this section, we establish a bound between the user and the system optimiza- 
tion problems (i.e. decentralized versus centralized problems) in the context of 
minimizing cost. For the (UO) decentralized problem we will consider the vari- 
ational inequality problem of finding x u £ K satisfying 

F( x u ) t (x — x u ) > 0, for all x £ K. 

Let x u and x s denote solutions of the user and system optimization problems 
respectively. Let Z u = F{x u ) t x u be the total cost for the user-optimized problem 
(UO) and Z s = F(x s ) t x a = min x£ K F(xYx be the total cost for the system- 
optimized problem (SO). In this section, we provide a bound on Z u jZ s for cost 
functions F(x) = Gx + 6, with G >~ 0 (i.e. positive definite) and asymmetric 
matrix, tfx > 0 for all x £ K (notice that this follows when constant vector 
b > 0 and K C J?” ). In this case, the system optimization problem involves the 
minimization of a strictly convex quadratic function over the set K. 
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2.1 A Measure of Asymmetry of a Matrix 

For a matrix G, we consider the symmetrized matrix 

C. G + G t 
2 

and introduce the following measure c 2 of the degree of asymmetry of matrix G : 



Definition 1. 



c 



2 



\\s 1 G < ||5 = sup 

w^O 



l|S- 1 G W ||| 

Mil 



sup 

0 



w t G t S~ 1 Gw 

w t Sw 



\\S~ 1/2 GS~ 1/2 \\ 2 . 



Note that by setting l = S 1 / 2 w, the previous definition of c 2 becomes 

ItS-W&S-'GS-Wl 



c = sup ■ 
1^0 



= A max {S- 1 ' 2 G t S- l GS- 1 / 2 ). 



When the matrix G is positive definite and symmetric, that is, G = G* (and 
therefore, S = G), then c 2 = 1. As an example, consider 



G = 



1 a 
— a 1 



Since S = I, it easily follows that c 2 = 1 + o 2 . The quantity c 2 in this case 
quantifies the observation that as |a| increases, the degree of asymmetry of G 
increases as well. 

The next lemma indicates that if G 2 >: 0, then the degree of asymmetry of G is 
bounded. 



Lemma 1. (see [11]) If G 2 is a positive semidefinite matrix then c 2 = 
||MG|||<2. 



2.2 The Bound 

Theorem 1. For an affine variational inequality problem with problem function 
F(x) = Gx + b, with G >~ 0, b*x > 0 for all a : £ K, we have: 





if c 2 < 2 
if c 2 > 2. 



Proof: 

We first observe that since x u solves VI(F,K) and x s £ K then the variational 
inequality formulation implies that 
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Z u F^Xu) X u — F{^Xu) X s 

= x t u G t x s + tfxs 

= x t u G t S~ 1 Sx s + tfxs (multiplying with S and S' -1 ) 

< ||(x„) t G t 5 _1 ||5 ||x s ||s + tfxg (from Cauchy's inequality, and 5 >- 0) 

< H^mIIs ||-S' _1 G||s ||x g ||s + tfxg (from the norm inequality) 

= C ||x u ||s ||x s ||s + 6 4 x s . (Definition 1) 

For every aq, 02 > 0, we have 

0 < (VoT||a; u ||s - V^ 2 \\ x s\\s) 2 = aq||x„||| + a 2 ||x g ||| - 2 v / alai'||x u ||s||x s ||s 
leading to 

2 v /aio 2 ||x„||s||x s ||s < ai||x„||! + a 2 ||x g |||. 

This in turn implies that if we choose 2^/cq a 2 > c, then 

c||x u ||s||x s ||s < ai||x„||| + a 2 ||x s |||. 

2 

We thus obtain that for all aq, a 2 > 0 and aia 2 > 

Z u < ai||x u ||| + a 2 ||x g ||| + 6*x s 
= aix^Sxu + a 2 xlSx s + 6*x s 

= ai(x t u G t x u + tfx-u) + a 2 (x t s G t x s + 6*x s ) - aitfxu - (a 2 - l)6 4 x s . 

If we further select a 2 > 1 and since 6‘x u , 6*x s > 0, we obtain 

Z u "G oq Z u + a 2 Z s . 



If we further impose the condition aq < 1, we obtain the bound: 

Z u < a 2 

Z s ~ 1 - 01 ‘ 

Given that we have freedom to select aq and a 2 , we find the best upper bound 
by solving 

a 2 

minimize 

1 — aq 

subject to aia 2 > ^- ' > 

a 2 > 1 , 0 < ai < 1 . 

The optimal solution to Problem (1) is given as follows: 

If c 2 < 2, the optimal solution is aq = c 2 /4, a 2 = 1 leading to 

4 

Z a ~ 4 - c 2 ' 

If c 2 > 2, the optimal solution is aq = 1/2, a 2 = c 2 /2 leading to 




.□ 



Remark: Notice that the results in this section apply to a feasible region that 
is non-convex and as a result includes integers. 
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(a) Separable affine functions: 

When the variational inequality problem function F is separable, it has 
components Fi(x) = giXi + bi. In this case the matrix G is diagonal, with 
diagonal elements (ji >0. In this case c 2 = 1 and the bound in Theorem 1 
becomes 

Z u 4 4 

Z s - A- c 2 3 

originally obtained in Roughgarden and Tardos [21]. 

(b) Non-separable symmetric affine functions: 

When the variational inequality problem function F is non-separable, that 
is F(x) = Gx+b, with G a general symmetric positive definite matrix, then 
c 2 = 1 and thus Z u /Z s < 4/3. This illustrates that the bound of 4/3 by 
Chau and Sim [2] is also a special case of the bound in this paper. 

(c) Non-separable asymmetric affine functions: 

When the matrix G is “not too asymmetric” (in the sense that c 2 < 2) 
then the bound becomes 4 f c2 . On the other hand, for “rather asymmetric” 
matrices (in the sense that c 2 > 2 ) then the bound becomes c 2 . 



Corollary 1. (see [19]) When the constant term is zero in the variational in- 
equality problem function, that is, F( x) = Gx, where G is a positive definite 
asymmetric matrix, then 




We now discuss some examples that illustrate the tightness of the bounds 
we established. In this subsection we discuss some examples that illustrate the 
tightness of the bounds we established. In particular through the example below 
we illustrate how the bound becomes tight for some families of problems. 

Example: 

Consider F\(x) = ax 1 + fx 2 + b±, F 2 (x) = — fx 1 + gx 2 + 62 and feasible region 
K = {x = (x 1 , x 2 ) : x 1 + x 2 = d, x 1 , x 2 > 0}. 

Notice that the degree of asymmetry constant c 2 = A maa; (S' _ 1 / 2 G t S' _ 1 GS' _1 / 2 ) = 
1 + £, si » cetathisH!m , pleG= (_“ / /), s = 

We will consider the following cases. 

(a) First we consider the case where a < g. Let / = a and constant terms 
b\ = b = (g — a)d and b 2 = 0. Notice that the measure of asymmetry 
becomes c 2 = 1 + - < 2, since - < 1. Furthermore, since a < q then 
bi = b > 0. The user-optimized solution is X* = 0, x\ = d and Z u = gd 2 . 
The system-optimized solution is x\ = 2 (a+g) — 0 (since b = (g — a)d), 

x 2 a = 2 (a+g) — 0 an d Z s = _ Notice that then 

Z u ^g 4 4 

Z s 3 g — a 3 — | 4 — c 2 
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Furthermore, when in particular, a = g then c 2 = 2 and ^ ^ „ = 2 = 

c 2 . 

(b) Consider the case where a > g. Choose as / = —a and as constant terms 

b\ = 0, b 2 = b = (a — g)d. Notice that since a > 5, 62 = b > 0. Moreover, 
the measure of asymmetry becomes c 2 = 1 + - > 2. The user-optimized 
solution is x* = d, x 2 = 0 and Z u = ad 2 . The system-optimized solution 
is x l = > 0- x2 s = IZaTg) - 0 ( since b 2 = b = (a - g)d) and 

Z s = a ( 2(a+g ) ) 2 + 9( 2 (a+g) ) 2 + b 2(a+g) ■ Suppose that d -> oo. Notice that 
then as d — » oo, 

-Zu ^ a + g _ 1 a _ c 2 

Z s g g 

These two examples establish that the bound in Theorem 1 is tight. Fur- 
thermore, this last example also suggests that the bound is tight even when 
the demand d is very large. 

(c) Finally, consider the case where the constant terms b\ = b 2 = 0. Choose as 
/ = —a. Then c 2 = 1 + |. The user-optimized solution is x\ = d, a: 2 = 0 and 

Z u = ad 2 . The system-optimized solution is x\ = ^+ g ) > 0, x 2 = ^+ g ) — 0 

and Z s = 23^-' Notice that then 

* a+g 



Z u a + g a 2 

— = = 1 + - = c . 

Z s g g 

This example establishes that the bound in Corollary 1 is tight (i.e., when 
the constant term is zero). 



In this example G = 



a f 
-f 9 



metrized version is matrix S = G+G = 



is a positive definite matrix since its sym- 
Furthermore, matrix 



a 0 
0 9 



G 2 becomes ^ ^ ^ ^ a J 2 ^ with a symmetrized matrix G + j, G ^ = 

a 2 - f 2 0 \ 

q f' 2 )' ^ nce * n f am ily °f examples above, we chose f = a 

or / = —a, we observe that in both cases the symmetrized matrix of G 2 is 



0 0 
0 g 2 — a 2 



Therefore, G 2 is a positive semidefinite matrix when a < g. This 



is equivalent to w t G 2 w = (G t w) t (Gw) > 0, which suggests that the angle be- 
tween vectors ( G t w ) and (Gw) is less than or equal to 90 degrees. Notice that 

< 2. Furthermore, G 2 is a negative 



in this case it is also true that c 2 



= 1 + f 



definite matrix when a > g. Since w t G 2 w = (G t w) t (Gw) < 0, it follows that 
the angle between vectors (G t w) and (Gw) is more than 90 degrees. It is also 
the case that c 2 = 1 -I- | >2. This latter observation illustrates further how the 
bound on constant c 2 connects with the degree of asymmetry of matrix G. 
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3 A Bound for Nonlinear, Asymmetric Functions 

In this section, we assume that the Jacobian matrix is not a constant matrix 
G but a positive definite, nonlinear and asymmetric matrix \7F(x). The posi- 
tive definiteness assumption of the Jacobian matrix implies that the variational 
inequality problem has a unique solution (see for example [17] for details). 

We introduce the symmetrized matrix, S(x ) = ^f(x)+vf{x) . now extend 
Definition 1 for measuring the degree of asymmetry of the problem to a 
definition that also incorporates the nonlinearity of the cost functions involved. 

Definition 2. We define a quantity c 2 that measures the degree of asymmetry 
of the Jacobian matrix \ / F(x). That is, 

C 2 = sup ||5(a;)- 1 VF(a;)|| 2 (:c) . 

xGK 

Constant c 2 is in this case the supremum over the feasible region, of the maximum 
eigenvalue of the positive definite and symmetric matrix 

S{x)- l/2 \/F{x) t S{x)- 1 \/F{x)S{x)- 1/2 . 

When the Jacobian matrix is positive definite and symmetric, then c 2 = 1. 
Furthermore, we need to define a measure of the nonlinearity of the problem 
function F. As a result, we consider a property of the Jacobian matrix which 
always applies to positive definite matrices. This allows us in some cases to 
provide a tight bound. This bound naturally extends the bound in Theorem 
1 from affine to nonlinear problems. The bound involves the constant A that 
measure the nonlinearity of the problem. 

Definition 3. (see [23] for more details ) 

The variational inequality problem function F : R n — > R n satisfies Jacobian 
similarity property if it has a positive semidefinite Jacobian matrix (VF(x) V 
0, \/x £ I\ ) and Vie € R n , and Mx, x € K, there exists A > 1 satisfying 

-j w t VF(x)w < w t X7F(x)w < Aw t X7F(x)w. 

l 



Lemma 2. (see [23]) The Jacobian similarity property holds under either of the 
following conditions: 

— The Jacobian matrix is strongly positive definite (i.e. has eigenvalues 
bounded away from zero). Then a possible bound for the constant A is 

^ _ max xeK A max (S(x)) 

A min (S(x)) ' 

— The problem function is affine, with positive semidefinite Jacobian matrix 
G. In this case A = 1. 
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Theorem 2. (see [19]) For a variational inequality problem with a strictly 
monotone, nonlinear continuously differentiable problem function F satisfying 
the Jacobian similarity property, F(0Y x > 0 for all x £ I\ , we have: 




4 

4 — c 2 A 

c 2 A 2 - 2 (A - 1) 



*/c 2 <! 

*/ ° 2 > A - 



Remarks: 



1. When the variational inequality problem function is affine, F(x) = Gx + b, 
then X/F(x) = G and as a result Definition 3 holds for A = 1. Therefore, 
= 2 and the bound coincides with the one we found in the previous section. 



2. When the variational inequality problem function has a symmetric Jacobian 
matrix, as we discussed c 2 = 1. Therefore, 

i. If A < 2 it follows that c 2 = 1 < and < jYj- 

ii. If A > 2 it follows that c 2 = 1 > and ^ < 1 + A 2 . 

3. When the term E(0) = 0 then similarly to Corollary 1, the bound becomes 

j^<c 2 A 2 ~2(A~l). 

We discuss an example where the bound we established for nonlinear problems 
is tight. 

Example: 

Consider the feasible region 

K = {x = (a; 1 , x 2 ) £ R\ : x 1 + x 2 = 1/2}. 

Consider cost functions 



Fi(x) = {x 1 ) 2 + x 1 + ax 2 , Fffx) = {x 2 ) 2 + x 2 — ax 1 . 



The Jacobian matrix is 



VF(x) 



( 2a: 1 + 1 a \ 

^ —a 2x 2 + 1 J ' 



In this example A = maXxeJf = 1 = 2. The constant c 2 = 1 + a 2 . The 

system-optimized solution is x s = (1/2, 1/2) which gives rise to Z s = F(x s Yx s = 
|. Furthermore, the user-optimized solution is a: 1 = — p, a; 2 = -j 2 . Then 
Z u = F(x u Yx u = ^ As a result, 

5a 3 5a 2 a 

= ~^4 + U ~ 2 + ' 

The bound in this case becomes -p- < c 2 A 2 — 2 (A — 1) = 4(1 + a 2 ) — 2 = 2 + 4a 2 . 
From this discussion we conclude that the bound is tight for all values of a 
satisfying 

5a 3 5a 2 a 0 q 0 

— — — + — — — + 1 = 2 + 4a 2 <+> 5a 3 + 86a 2 + 12a + 24 = 0. 

24 12 2 
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4 Bounds via Semidefinite Optimization 



In this section, we illustrate how to bound the “price of anarchy” using ideas from 
semidefinite optimization. We illustrate a method to derive the same bounds as 
in the previous sections when the Jacobian matrix of the problem function, is 
not necessarily invertible, i.e., it is a positive semidefinite rather than a positive 
definite matrix. 



Theorem 3. For an affine variational inequality problem with problem function 
F(x) = Gx + b, with G F 0, b*x > 0 for all x € K, the optimal objective function 
value of the following semidefinite optimization problem, 



a 2 

mm 

1 — Oi 



(2) 



satisfying 




>: 0 



a 2 > 1, 0 < oi < 1, 
determines the bound < , 02 . . 

Z s — 1 — a? 



Proof: 

As in Theorem 1 since x u solves the variational inequality formulation and x s £ 
K then 

Ffxufxu < Ffxufxs = (. x u ) t G t x s + tfxs. (3) 

We wish to find oi, a 2 > 0, so that 

x t u G t x s < aix^Sxu + a 2 xlSx s (since S = G+ 2 G ) 

= aix^Gxu + 02 xlGx s 



Notice that (4) is equivalent to 








> 0 . 



This relation follows when the symmetric part of the 2 n x 2 n matrix 

( a 'fa S') - 0 (Ulat (“4 i ) - 0) ' 

Furthermore, if 02 > 1 and oi > 0 then since b*x u , tfxg > 0, relations (3), (4) 
imply that 

Z u — F'fxy,') A Oi F{x u f x u T a 2 F(x s ) x s . 



Z u ^ a 2 

Z s ~ 1 - ai 



Therefore, if a\ < 1 then 
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Given that we have freedom to select ai and 02 , we find the best upper bound 
by solving the following semidefinite optimization problem, 

° 2 

mm (5) 

1 — a\ 



satisfying 




>: 0 . 



a 2 > 1, 0 < or < 1. 



Notice that the optimal objective function value to the semidefinite optimization 
problem (5) is the bound determining the “price of anarchy”. □ 

Remark: This proof does not require the matrix G to be positive definite. 
Matrix G in this case is positive semidefinite. The next corollary connects the 
bounds in Theorem 3 and 1 illustrating that when the matrix G 0, the two 
bounds coincide. 



Corollary 2. For an affine variational inequality problem F( x) = Gx + b, with 
matrix G >~ 0, Theorems 1 and 3 yield the same bound. 



Proof: 

Notice that when the matrix G 0 (and as a result, its symmetrized matrix 



S >~ 0) then the matrix 




^ 0 if and only if ai > 0 and the matrix 



( ail 

! 1 1 

\ -S' 2 GS~ 2 
\ 2 



1 . 1 

-S" 2 G t S~2 
2 



a 2 X 



>= 0 , 



where matrix X is the identity matrix. Notice that this is also in turn equivalent 

1 . , 1 
rr — ^ /-it n — 1 /-( n — 

to ai > 0 and the matrix a\ail — - — - — — - being positive semidefinite. 

x (s~bG t s- 1 GS- 1 



But this is equivalent to ai > 0 and 0102 > 



= x ( see 



4 “4 

Definition 1). As a result the SDP involved in Theorem 3 yields in this case the 

same bound as the optimization problem (1) in Theorem 1. □ 

Remark: Notice that when the constant term is zero then the SDP problem 
reduces to 

°2 , a \ 

mm ( 0 ) 

1 — CL 1 

satisfying ^ j ^0 



a 2 > 0, 0 < ai < 1. 

Similarly to Corollary 2, one can show that when G >~ 0, the bound is c 2 . Notice 
that this is the same bound as in Corollary 1. 

This approach also extends to nonlinear problems. 
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Theorem 4. (see [19]) For a variational inequality problem with a positive 
semidefinite Jacobian matrix (i.e., \7F(x) F OJ satisfying the Jacobian simi- 
larity property, F(0)*x > 0, for all x € K, the optimal objective function value 
of the following semidefinite optimization problem, 



a 2 

mm 

1 — ai 



( 7 ) 



satisfying „ X) ) “ °' ^ 6 * 

a 2 > 1, 0 < ai < 1, 

determines the bound 4^ < , " 2 . . 

Z s — 1— a J 

Corollary 2 extends to nonlinear problems as well. 



Corollary 3. (see [19]) When the Jacobian matrix is strongly positive definite 
then the bound in Theorem f becomes the same as in Theorem 2. That is, 




4 

4 — c 2 A 

c 2 A 2 - 2 (A - 1) 



*/c 2 <| 
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Abstract. We study the problem of minimizing the maximum latency of 
flows in networks with congestion. We show that this problem is NP-hard, 
even when all arc latency functions are linear and there is a single source 
and sink. Still, one can prove that an optimal flow and an equilibrium flow 
share a desirable property in this situation: all flow-carrying paths have 
the same length; i.e., these solutions are “fair,” which is in general not 
true for the optimal flow in networks with nonlinear latency functions. 
In addition, the maximum latency of the Nash equilibrium, which can 
be computed efficiently, is within a constant factor of that of an optimal 
solution. That is, the so-called price of anarchy is bounded. In contrast, 
we present a family of instances that shows that the price of anarchy is 
unbounded for instances with multiple sources and a single sink, even in 
networks with linear latencies. Finally, we show that an s-f-flow that is 
optimal with respect to the average latency objective is near optimal for 
the maximum latency objective, and it is close to being fair. Conversely, 
the average latency of a flow minimizing the maximum latency is also 
within a constant factor of that of a flow minimizing the average latency. 



1 Introduction 

We study static network flow problems in which each arc possesses a latency 
function, which describes the common delay experienced by all flow on the arc 
as a function of the flow rate. Load-dependent arc costs have a variety of ap- 
plications in situations in which one wants to model congestion effects, which 
are bound to appear, e.g., in communication networks, road traffic, or evacua- 
tion problems. In this context, a unit of flow frequently denotes a huge number 
of “users” (or “agents”), which might represent data packages in the Internet, 
drivers on a highway system, or individuals fleeing from a disaster area. Depend- 
ing on the concrete circumstances, the operators of these networks can pursue a 
variety of system objectives. For instance, they might elect to minimize the aver- 
age latency, they might aim at minimizing the maximum latency, or they might 
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try to ensure that users between the same origin-destination pair experience es- 
sentially the same latency. In fact, the ideal solution would be simultaneously 
optimal or near optimal with respect to all three objectives. 

For linear latencies, we prove the existence of an s-f-flow that is at the same 
time optimal for two of the three objectives while its average latency is within a 
factor of 4/3 of that of an optimum. As attractive as this solution might be, we 
also show that it is NP-hard to compute. Moreover, there is a surprising differ- 
ence between linear and nonlinear latency functions. Namely, this particular flow 
remains optimal with respect to the maximum latency and near optimal with 
respect to the average latency, but it does in general not guarantee that different 
users face the same latency. However, an optimal s-t-flow for the average latency 
objective can be computed in polynomial time, and we show that the latency 
of any one user is within a constant factor of that of any other user. In partic- 
ular, the maximum latency is within the same constant factor of the maximum 
latency of an optimal solution to the latter objective. This constant factor only 
depends on the class of allowable latency functions. For instance, its value is 2 
for the case of linear latencies. Linear latencies are sufficient for certain conges- 
tion phenomena to occur. One interesting example is Braess’ paradox (1968), 
which refers to the fact that the addition of an arc can actually increase the 
(average and maximum) latency in a network in which users act selfishly and 
independently. This user behavior is captured by the Nash equilibrium of the 
underlying game in which each user picks a minimal latency path, given the 
network congestion due to other users (Wardrop 1952). While the inefficiency 
of this so-called user equilibrium and hence the severity of Braess’ paradox had 
previously been bounded in terms of the average latency, it turns out that it 
is also bounded with respect to the maximum latency. Indeed, the latencies en- 
countered by different users between the same origin-destination pair are the 
same. The user equilibrium therefore represents another flow that can be com- 
puted in polynomial time and that is optimal or close to optimal for the three 
objectives introduced earlier. 

The Model. We consider a directed graph G = (N,A) together with a set of 
source-sink pairs K C TV x N. For each terminal pair k = ( Sk,tk ) G K, let 
Vk be the set of directed (simple) paths in G from Sk to tk, and let dk > 0 
be the demand rate associated with commodity k. Let V := (J keK Pk be the 
set of all paths between terminal pairs, and let d := Yhk^K dk be the total 
demand. A feasible flow / assigns a nonnegative value fp to every path P £ V 
such that Ylppv k f p = dk for all k € K. In the context of single-source single- 
sink instances, we will drop the subindex k. Each arc a has a load-dependent 
latency £ a (-)- We assume that the functions l a : IR>o IR>o are nonnegative, 
nondecreasing, and differentiable. We define the latency of a path P £ V under 
a given flow / as £p(f) := YlaeP ^(Y^QeV:Q 3 a /q)- The maximum latency of 
a feasible flow / is L(f ) := max{tp(/) : P £ V , fp > 0}. We call a feasible 
flow that minimizes the maximum latency a min-max flow and denote it by /. 
The maximum latency problem consists of finding a min-max flow. The average 
latency of a feasible flow / is defined as C{f) := Ylp^v ^p(f)fp/d- We refer 
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to the optimal solution with respect to this objective function as the system 
optimum and denote it by /*. A feasible flow is at Nash equilibrium (or is a 
user equilibrium ) if for every k £ K and every two paths Pi,P 2 £ Vk with 
fp 1 > 0, £p 1 (f) < fp 2 (/). In other words, all flow-carrying Sk~tk ~ paths have 
equal (and actually minimal) latency. In particular, equilibrium flows are “fair,” 
i.e., they have unfairness 1, if the unfairness of a feasible flow / is defined as 
max fc6A -max{f Pl (/)/£p 2 (/) : P 1 ,P 2 £ Vk,fp 1 Jp 2 > 0}. 

Main Results. While a user equilibrium can be computed in polynomial time 
(Beckmann, McGuire, and Winsten 1956), and so can a system optimum if 
x£ a {x) is convex for all arcs a £ A, we show in Section 2 that it is an NP- 
lrard problem to compute a min-max flow. This result still holds if all latencies 
are linear and there is a single source-sink pair. Note that the flows that we 
are considering are not required to be integer, neither on paths nor on arcs. As 
pointed out earlier, a Nash equilibrium has unfairness 1 by construction. In Sec- 
tion 3, we establish the somewhat surprising existence of a min-max flow that 
is fair too, when latencies are linear and there is a single source and a single 
sink. In addition, although it is well known that system optima are unfair, we 
provide a tight bound that quantifies the severity of this effect. This bound ap- 
plies to general multicommodity flows and arbitrary latency functions. Finally, 
in Section 4, we show that in the single-source single-sink case under arbitrary 
latency functions, there actually exist solutions that are simultaneously optimal 
or near optimal with respect to all three criteria (maximum latency, average 
latency, and unfairness). In fact, this property is shared by the system optimum, 
the user equilibrium, and — to some extent — the min-max flow, albeit with dif- 
ferent bounds. Table 1 presents the bounds obtained for the three criteria in the 
single-source single-sink case with linear latencies. An important consequence of 
these results is that computing a user equilibrium or a system optimum con- 
stitutes a constant-factor approximation algorithm for the NP-hard maximum 
latency problem. On the other hand, already in networks with multiple sources 
and a single sink, the ratio of the maximum latency of a Nash equilibrium to 
that of the min-max flow is not bounded by a constant, even with linear latency 
functions. 

Related Work. Most papers on evacuation problems consider constant 
travel times; we refer the reader to the surveys by Aronson (1989) and 
Powell, Jaillet, and Odoni (1995) for more details. One exception is the 
work by Kohler and Skutella (2002). They considered a dynamic quickest 
flow problem with load-dependent transit times, for which they estab- 
lished strong NP-hardness. They also provided an approximation algorithm 
by considering the average of a flow over time, which is a static flow. 
Kohler, Langkau, and Skutella (2002) proposed to use time-expanded networks 
to derive approximation algorithms for a similar problem. 

The concept of the price of anarchy, which is the ratio of the per- 
formance of a Nash equilibrium to that of an optimal solution, was 
introduced by Koutsoupias and Papadimitriou (1999) in the context of 
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Table 1 . Summary of results for single-source single-sink networks with linear 
latency functions. The first entry in each cell represents a worst-case bound on 
the ratio of the value of the flow associated with the corresponding row to the 
value of an optimal flow for the objective function denoted by the corresponding 
column. The second entry refers to the theorem in this paper in which the 
respective result is proved. All bounds are tight, as examples provided after each 
theorem demonstrate. The bound of 4/3 on the ratio of the average latency of the 
user equilibrium to that of the system optimum was first proved by Roughgarden 
and Tardos (2002); we give a simpler proof in Theorem 7. Weitz (2001) observed 
first that this bound carries forward to the maximum latency objective for the 
case of only one source and sink; we present a generalization of this observation 
to multicommodity flows in Theorem 8. 





maximum latency 


average latency 


unfairness 


min-max flow 


1 


4/3 Thm. 11 


1 Thm. 5 


system optimum 


2 Thm. 10 


1 


2 Thm. 6 


Nash equilibrium 


4/3 Thm. 8 


4/3 Thm. 7 


1 



a game motivated by telecommunication networks. This inspired con- 
siderable subsequent work, including Mavronicolas and Spirakis 2001; 
Koutsoupias, Mavronicolas, and Spirakis 2003; Czumaj and Vocking 2002; 
Czumaj, Krysta, and Vocking 2002. These papers study the maximum latency 
of transmissions in two-node networks consisting of multiple links connecting 
a single source with a single sink. Under certain assumptions, the maximum 
latency when users are selfish is not too large compared to the best coordinated 
solution. Although these results are similar in nature to some of ours, the two 
models are not comparable because they work with a finite number of players 
and consider mixed strategies. In contrast, in our setting, every player just 
controls an infinitesimal amount of flow, making mixed strategies essentially 
irrelevant. Moreover, we work with arbitrary networks. For more details on the 
various routing games, we refer the reader to the survey by Czumaj (2004). 

Roughgarden and Tardos (2002), Roughgarden (2003), Schulz and Stier 
Moses (2003), and Correa, Schulz, and Stier Moses (2003) studied the price of 
anarchy with respect to the average travel time in general networks and for 
different classes of latency functions. In particular, if C is the set of allowable 
latency functions, the ratio of the average travel time of a user equilibrium to 
that of a system optimum is bounded by a(£), where a(£) is a constant that 
only depends on C. For example, in the case that C only contains concave func- 
tions, a(C) = 4/3. We will later make use of this result (Section 4). For the 
maximum latency objective, Weitz (2001) was the first to observe that the price 
of anarchy is bounded in single-source single-sink networks. He also presented a 
family of examples that showed that Nash equilibria can be arbitrarily bad in 
multiple commodity networks. Roughgarden (2004) gave a tight bound for the 
single-source single-sink case that depends on the size of the network. 





Computational Complexity, Fairness, and the Price of Anarchy 



63 



Game-theoretic concepts offer an attractive way of computing approximate 
solutions to certain hard problems. Indeed, Schulz and Stier Moses (2003) 
computed a provably good Nash equilibrium in a setting with mul- 
tiple equilibria in which computing the best equilibrium is hard. 
Anshelevich, Desgupta, Tardos, and Wexler (2003) approximated optimal solu- 
tions to a network design problem that is NP-lrard with the help of Nash and 
approximate Nash equilibria. A related idea was used by Fotakis et al. (2002) 
and Feldmann et al. (2003) to show that although it is hard to find the best and 
worst equilibrium of the telecommunication game described before, there exists 
an approximation algorithm for computing a Nash equilibrium with minimal 
social cost. 

In the context of Section 3, we should point out that there exist mul- 
tiple (nonequivalent) definitions of (un)fairness. The definition we use here 
comes from the competition between different agents in the routing game. 
Roughgarden (2002) introduced a less pessimistic version of unfairness, namely 
the ratio of the maximum latency of a system optimum to the latency of a user 
equilibrium; we later recover the bound that he obtained as a corollary to a 
more general result. Jahn, Mohring, Schulz, and Stier Moses (2002) considered 
the definition of unfairness presented here; they looked at flows that minimize 
the total travel time among those with bounded unfairness. 

2 Computational Complexity 

In our model, both the system optimum and the Nash equilibrium can 
be computed efficiently because they represent optimal solutions to cer- 
tain convex programs. 1 On the other hand, it follows from the work 
of Kohler and Skutella (2002) on the quickest s-f-flow problem with load- 
dependent transit times that the maximum latency problem considered here 
is NP-hard (though not necessarily in NP) when latencies include arbitrary non- 
linear functions or when there are explicit arc capacities. Lemma 1 below implies 
that the general maximum latency problem is in NP, while Theorem 3 estab- 
lishes its NP-lrardness, even in the case of linear latencies and a single source 
and a single sink. Note that the following result does not follow from ordinary 
flow decomposition as it is not clear how to convert a flow on arcs into a path 
flow such that the latency of the resulting paths remains bounded; in fact, it is 
a consequence of Theorem 3 that the latter problem is NP-hard, too. 

Lemma 1. Let f be a feasible flow for a multicommodity flow network with 
load- dependent arc latencies. Then there exists another feasible flow f such that 
L(f') < L(f), and f uses at most |A| paths for each source-sink pair. 

Proof. Consider an arbitrary commodity k € K. Let Pi, . . . ,P r be Sk-tk- paths 
such that fp i > 0 for i = 1 ,...,r, and /p 4 = ( h- Slightly overloading 

notation, we let Pi, ... , P r also denote the arc incidence vectors of these paths. 

1 This statement is only true for the system optimum if we assume that C(f) is convex. 
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Let’s assume that r > |A|. (Otherwise we are done.) Hence, the vectors P\,. . . ,P r 
are linearly dependent and ^ ?? = 0 h as a nonzero solution. Let’s assume 

without loss of generality that A r yf 0. We define a new flow f" (not necessarily 
feasible) by setting fp. := fp i — fp r for * = 1, ... , r, and fp := fp for all other 
paths P. Notice that under /", the flow on arcs does not change: 

r r — 1 i — 1 . r 

E p ^'k = E - E y p ifPr = E ■ 

i = 1 i—1 i=l Ar i = 1 

Here, we used the linear dependency for the last equality. In particular, L{f") < 
L(f). Let us consider a convex combination /' of / and f" that is nonnegative 
and uses fewer paths than /. Note that such a flow always exists because fp =0, 
and the flow on some other paths P±, . . . , P r -i might be negative. Moreover, 
L(f') < L(f), too. If /' still uses more than \A\ paths between Sk and tk, we 
can iterate this process so long as necessary to prove the claim. □ 



Corollary 2. The decision version of the maximum latency problem is in NP. 

Proof. Lemma 1 shows the existence of a succinct certificate. Indeed, there is a 
min- max flow using no more than \K\ ■ |A| paths. □ 

We are now ready to prove that the maximum latency problem is in fact 
NP-hard. We present a reduction from Partition: 

Given: A set of n positive integer numbers q\,...,q n . 

Question: Is there a subset / C {1, . . . , n} such that J2 ieI ft = ft? 

Theorem 3. The decision version of the maximum latency problem is NP- 
complete, even when all latencies are linear functions and the network has a 
single source-sink pair. 

Proof. Given an instance of Partition, we define an instance of the maximum 
latency problem as follows. The network consists of nodes 0,1,..., n with 0 
representing the source and n the sink. The demand is one. For i = 1 ,... ,n, 
the nodes i — 1 and i are connected with two arcs, namely a, with latency 
t ai (x) = qix and a, with latency = qt- 

Let L := | q%- Notice that the system optimum /* has cost equal to L 

and f* = 1/2 for all a £ A. We claim that the given instance of Partition is a 
Y ES-instance if and only if there is a solution to the maximum latency problem 
of maximum latency equal to L. Indeed, if there is a partition I, the flow that 
routes half a unit of flow along the 0-n-patlr composed of arcs a*, i € /, and cq, 
i £ I, and the other half along the complementary path has maximum latency L. 

To prove the other direction, assume that we have a flow / of maximum 
latency equal to L. Therefore, C(f) < L (there is unit demand), which implies 
that C(f) = L (it cannot be better than the optimal solution). As the arc flows of 
a system optimum are unique, this implies that f a = 1/2 for all a £ A. Take any 
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path P such that fp> 0, and partition its arcs such that I contains the indices 
of the arcs a; G P. Then, f X )" =1 qi = L = £ P (f) = J2 ie i f + and 

subtracting the left-hand side from the right-hand side yields ^ 

□ 



Corollary 4. Let f be a (path) flow in an s-t-network with linear latencies. Let 
(f a ■ a G A) be the associated flow on arcs. Given just ( f a '■ a G A) and L(f), 
it is NP-hard to compute a decomposition of this arc flow into a (path) flow f 
such that L(f') < L(f). In particular, it is NP-hard to recover a min-max flow 
even though its arc values are given. 

Note that Corollary 4 neither holds for the system optimum nor the user 
equilibrium. Any flow derived from an ordinary flow decomposition is indeed 
an optimal flow or an equilibrium flow, respectively. Let us also mention that 
Theorem 4.3 in Kohler and Skutella (2002) implies that the maximum latency 
problem is APX-lrard when latencies can be arbitrary nonlinear functions. 

3 Fairness 

User equilibria are fair by definition. Indeed, all flow-carrying paths between the 
same source and sink have equal latency. The next result establishes the same 
property for min-max s-f-flows in the case of linear latencies. Namely, a fair min- 
max flow always exists. Therefore, the difference between a Nash equilibrium 
and a min-max flow is that the latter may leave paths unused that are shorter 
than the ones carrying flow, a situation that cannot happen in equilibrium. The 
following result is not true for nonlinear latencies, as we shall see later. 

Theorem 5. Every instance of the single-source single- sink maximum latency 
problem with linear latency functions has an optimal solution that is fair. 

Proof. Consider an instance with demand d and latency functions £ a (fa) = 
q a f a + v a , for a £ A. Among all min-max flows, let / be the one that uses 
the smallest number of paths. Let Pi, P 2 , . . . , Pk be these paths. Consider the 
following linear program: 



min z 






(1) 


s.t. J2 (M 


+ r °) - z for * = 1 >- 


• -,k, 


(2) 


aePi 


Ph3a 






k 

i=l 


- d 




(3) 


f Pi > o 


for i = 1, . 


..,k. 


(4) 



Note that this linear program has k+ 1 variables. Furthermore, by construction, 
it has a feasible solution with 2 = L(f), and there is no solution with 2 < L(f). 
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Therefore, an optimal basic feasible solution gives a min- max flow that satisfies 
with equality k of the inequalities (2) and (4). As fp i > 0 for all i because of 
the minimality assumption, all inequalities (2) have to be tight. □ 

A byproduct of this proof is that an arbitrary flow can be transformed into 
a fair one without increasing its maximum latency. In fact, just solve the corre- 
sponding linear program. An optimal basic feasible solution will either be fair or 
it will use fewer paths. In the latter case, eliminate all paths carrying zero flow 
and repeat until a fair solution is found. 

Notice that the min-max flow may not be fair for nonlinear functions. Indeed, 
the instance displayed in Figure 1 features high unfairness with latencies that are 
polynomials of degree p> 2. When a = (1 + e) p_1 and b = 2 — (^p) P — 5 for 

some e > 0 and 6 > 0 such that b > 1, the min-max flow routes units of flow 
along the “top-bottom” and “bottom-top” paths, respectively, and units of 
flow along the “top-top” path. It is not hard to see that this flow is optimal. 
Indeed, the “bottom-bottom” path is too long to carry any flow. Moreover, by 
symmetry, the “top-bottom” and “bottom-top” paths have to carry the same 
amount of flow. Therefore, the optimal solution can be computed by solving 
a one-dimensional minimization problem, whose only variable is the amount x 
of flow on the “top-top” path. The unique optimal solution to this problem is 
x = 577 ■ Let us compute the unfairness of this solution. The “top-top” path has 

latency equal to 2 (|jr§) P , which tends to (^) P 1 as e — ■> 0. The latency of the 
other two paths used by the optimum is equal to 2 — 5. Therefore, the unfairness 
of this min-max flow is arbitrarily close to 2 P . 




ax p + b ax p + b 



Fig. 1 . Fair min-max flows may not exist for instances with nonlinear latencies. 



A typical argument against using the system optimum in the design of route- 
guidance devices for traffic assignment is that, in general, it assigns some drivers 
to unacceptably long paths in order to use shorter paths for most other drivers; 
see, e.g., Beccaria and Bolelli (1992). The following theorem quantifies the sever- 
ity of this effect by characterizing the unfairness of the system optimum. It turns 
out that there is a relation to earlier work by Roughgarden (2002), who com- 
pared the maximum latency of a system optimum in a single-sink single-source 
network to the latency of a user equilibrium. He showed that for a given class 
of latency functions £, this ratio is bounded from above by 7 (£), which is de- 
fined to be the smallest value that satisfies £* a (x) < 7 (C)£ a (x) for all £ € C and 
all x > 0. Here, £*(x) := £ a (x ) + x£' a (x ) is the function that makes a system 
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optimum for the original instance a user equilibrium of an instance in which 
the latencies are replaced by £* (Beckmann, McGuire, and Winsten 1956). For 
instance, 7(polynomials of degree p) = p + 1. We prove that the unfairness of a 
system optimum is in fact bounded by the same constant, even for general in- 
stances with multiple commodities. The same result was independently obtained 
by Roughgarden (personal communication, October 2003). 

Theorem 6. Let f* be a system optimum in a multicommodity flow network 
with arc latency functions drawn from a class C. Then, the unfairness of f* is 
bounded from above by 7 (£). 

Proof. We will prove the result for the single-source single-sink case. The ex- 
tension to the general case is straightforward. As a system optimum is a user 
equilibrium with respect to latencies £*, there exists L* such that £p{f*) = L* 
for all paths P £ V with ff > 0. From the definitions of £* and we have 
that £ a (x) < t*(;r) < / y(£)£ a (x) for all x. Let Pi, Pi € V be two arbitrary paths 
with /p ,/p > 0. Hence, £p 1 (f*) < L* and £p 2 (f*) > L* /7(F). It follows that 

£pAn/tpln<i{L)- □ 

Notice that Theorem 6 implies Roughgarden’s earlier bound for the single- 
source single-sink case. Indeed, for a Nash equilibrium /, min{fp(/*) : P £ 
P,/p > 0} < min {£p(f) : P € V , fp > 0}. Otherwise, C(/*) > C(f), which 
contradicts the optimality of /*. In addition, the example shown in Figure 2 
proves that the bound given in Theorem 6 is tight. Indeed, it is easy to see that 
the system optimum routes half of the demand along each arc, implying that 
the unfairness is £*(d/2)/£(d/2). Taking the supremum of that ratio over d > 0 
and £ € C yields 7 (£). 



(■Ad/ 2 ) 




Fig. 2. Instance showing that Theorem 6 is tight. 



4 Price of Anarchy and Related Approximation Results 

Nash equilibria in general and user equilibria in particular are 
known to be inefficient, as evidenced by Braess’ paradox (1968). 
Koutsoupias and Papadimitriou (1999) suggested measuring this degrada- 
tion in performance, which results from the lack of central coordination, by 
the worst-case ratio of the value of an equilibrium to that of an optimum. 
This ratio has become known as the “price of anarchy,” a phrase coined by 
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Papadimitriou (2001). It is quite appealing (especially for evacuation situations) 
that in the routing game considered here, the price of anarchy is small; i.e., 
the selfishness of users actually drives the solution close to optimality. Recall 
that the user equilibrium results from everyone choosing a shortest path under 
the prevailing congestion conditions. Since a user equilibrium can be computed 
in polynomial time, this also leads to an approximation algorithm for the 
maximum latency problem. 

In order to derive a bound on the price of anarchy for the maxi- 
mum latency objective, we use a corresponding bound for the average la- 
tency of Nash equilibria, which was first proved for linear latency functions 
by Roughgarden and Tardos (2002), and then extended to different classes 
of latency functions by Roughgarden (2003), Schulz and Stier Moses (2003), 
and Correa, Schulz, and Stier Moses (2003). For the sake of completeness, 
let us include a simpler proof of Roughgarden and Tardos’ result (see also 
Correa, Schulz, and Stier Moses 2003). 

Theorem 7 (Roughgarden and Tardos 2002). Let f be a user equilibrium 
and let f* be a system optimum in a multicommodity flow network with linear 
latency functions. Then C(f) < | C(f*). 

Proof. Let h a {x) = q a x + r a with q a , r a > 0 for all a G A. W.l.o.g., d = 1. Then, 
C{f) = 53 ( q a fa + r a )fa < 53(«a/a + r a )f* < 

aGA aGA 

5 >»/«* + r a )f : + \ 53 qafa < c(n + \cu) ■ 

aGA aGA 

The first inequality holds since the equilibrium flow / uses shortest paths with 
respect to the arc latencies caused by itself. The second inequality follows from 

(/:-/»/ 2) 2 >o. ' ' □ 

For general latency functions 

C(f) < a(C)C(f*), where a(£) := (l- sup ’(5) 

^ <e£,o<i<ii ^ d£(d) >) 

and the proof is similar to the one given above for Theorem 7; see 
Roughgarden (2003) and Correa, Schulz, and Stier Moses (2003) for details. For 
polynomials with nonnegative coefficients of degree 2, a(C) equals 1.626; for 
those with degree 3, a{C) = 1.896; in general, a(C) = 0{jp/\np) for polynomials 
of degree p. It was first noted by Weitz (2001) that in networks with only one 
source and one sink, any upper bound on the price of anarchy for the average 
latency is an upper bound on the price of anarchy for the maximum latency. We 
include a multicommodity version of this result. 

Theorem 8. Consider a multicommodity flow network with latency functions 
in C. Let f be a Nash equilibrium and f a min-max flow. For each commodity 
k € K, Lk(f) < -^fa{C)L{f), where Lk is the maximum latency incurred by 
commodity k, dk is its demand rate, and d is the total demand. 
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Proof. Let f* be the system optimum. Then, 

d k L k (f) < dC(f) < da(C)C(f*) < da(C)C(f) < da(C)L(f) . 

Here, the first inequality holds because / is a Nash equilibrium, the second 
inequality is exactly (5), the third one comes from the optimality of /*, and the 
last one just says that the average latency is less than the maximum latency. □ 

The proof of Theorem 8 implies that if for a given single-source single-sink 
instance an equilibrium flow / happens to be a system optimum, then / is 
also optimal for the maximum latency objective. 2 For single-source single-sink 
networks with general latency functions, Theorem 8 shows that computing a 
Nash equilibrium is an a(£)-approximation algorithm for the maximum latency 
problem. Notice that this guarantee is tight as shown by the example given in 
Figure 3, which goes back to Braess (1968). Indeed, the latency of a Nash equi- 
librium is £(d) while the maximum latency of a min-max flow, which coincides 
with the system optimum, is 1(d) — ^max^{| (£(d) — £(x))}. Taking the supre- 

mum over d > 0 and £ € C, the ratio of the latency of the Nash equilibrium to 
that of the min-max flow is arbitrarily close to a(C). 




Fig. 3. Instance showing that Theorem 8 is tight for single-commodity networks. 



For instances with multiple sources and a single sink, the maximum latency 
of a user equilibrium is unbounded with respect to that of an optimal solution, 
even with linear latencies. In fact, we will show that the price of anarchy cannot 
be better than fi(n), where n is the number of nodes in the network. Note that 
this also implies that the price of anarchy is unbounded in single-source single- 
sink networks with explicit arc capacities. Weitz (2001) showed that the price of 
anarchy is unbounded in the case of two commodities, and Roughgarden (2004) 
proved that it is at most n — 1 if there is a common source and sink. 

Theorem 9. The price of anarchy in a single- commodity network with midtiple 
sources and a single sink is l?(n), even if all latencies are linear functions. 

Proof. Fix a constant £ > 0 and consider the instance presented in Figure 4. 
Nodes n, n— 1, . . . , 1 are the sources while node 0 is the sink. Nodes i and i — 1 are 

2 For instance, if all latency functions are monomials of the same degree but with 
arc-dependent coefficients, it is well known that user equilibria and system optima 
coincide (Dafermos and Sparrow 1969). 
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connected with two arcs: a, with constant latency equal to 1 and a* with latency 
equal to x/e % . Let the demand entering node i > 0 be £*. The user equilibrium 
of this instance routes the flow along paths of the form a,, Oj_i, . . . , a\ and has 
maximum latency n. To show the claim, it suffices to exhibit a good solution. 
For instance, for origin i, let its demand flow along the path a*, fij-i, . . . , «i. 
Under this flow, the load of cq is equal to £* +1 + •••+£" and its traversal time 
is (£ l+1 + • • • + £ n )/e l = e l + • • • + £ n ~ l . Hence, we can bound the maximum 
latency from above by 1 + which tends to 1 when £ — » 0. □ 
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Fig. 4. Nash equilibria can be arbitrarily bad when multiple sources are present. 

In the single-source single-sink case, Nash equilibria are not the only good 
approximations to the maximum latency problem; an immediate corollary of 
Theorem 6 is that system optima are also close to optimality with respect to the 
maximum latency objective. 

Theorem 10. For single-source single- sink instances with latency functions 
drawn from C, computing a system optimum is a 7 (£)- approximation algorithm 
for the maximum latency problem. 

Proof. Theorem 6 states that the length of a longest path used by the system 
optimum /* is at most 7 (£) times the length of a shortest flow-carrying path. 
The latter value cannot be bigger than the maximum latency of a path used 
by the min-max flow because /* is optimal for the average latency; the result 
follows. □ 

The bound given in Theorem 10 is best possible. To see this, consider the 
instance depicted in Figure 5. The min-max flow routes the entire demand along 
the lower arc, for a small enough £ > 0. On the other hand, the unique system 
optimum has to satisfy £*(f*) = £*(d) — e, where f* is the flow along the lower 
arc. Therefore, the upper arc has positive flow and the maximum latency is 
£*{d) — £. The ratio between the maximum latencies of the two solutions is 
arbitrarily close to £*(d)/£(d). Taking the supremum over d > 0 and £ € C 
shows that the bound in Theorem 10 is tight. 

To complete Table 1, let us prove that the average latency of the min-max 
flow is not too far from that of the system optimum. 

Theorem 11. Let f be a min-max flow and let f* be a system optimum for an 
instance with a single source, a single sink and latencies drawn from C. Then, 
C(f) < a(£)C(f*). 
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(*(d) — e 




Fig. 5. Instance showing that Theorem 10 is tight. 



Proof. Note that C(f) < L(f) < L(f) = C(f) < a(£)C(f *), where / is the 
Nash equilibrium of the instance. □ 

Again, the guarantee given in the previous theorem is tight. To show this, 
it is enough to note that the equilibrium flow and the min-max flow coincide in 
the example of Figure 6, and their average latency is ((d). Moreover, the average 
latency of the system optimum is arbitrary close to ((d) — | (( (d) — i(x )) } . 

Taking the supremum of the ratio of these two values over d > 0 and l £ £ 
completes the argument. 



((d) + e 




Fig. 6. Instance showing that Theorem 11 is tight. 



5 Concluding Remarks 

While Table 2 summarizes our findings for single-source single-sink networks with 
latencies drawn from a given class £ of allowable latency functions, let us briefly 
discuss another useful objective function. A different way of interpreting the 
two-node scenario of Koutsoupias and Papadimitriou (1999) in which each path 
consists of a single arc, is to minimize the maximum latency over flow-carrying 
arcs. In contrast to the maximum latency problem, this problem can be solved 
in polynomial time when latency functions are convex. However, the optimal 
“bottleneck flow” can be arbitrarily bad with respect to all three objectives 
considered in the main part of this paper. On the other hand, the equilibrium 
flow, the system-optimal flow, and the min-max flow can all at once be arbitrarily 
bad with respect to the bottleneck objective. We leave the details to the full 
version of this paper. 
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Table 2. Overview of approximation guarantees for single-source single-sink 
networks when latencies belong to a given set C. All bounds are tight. The “?” 
indicates that no upper bound is known; recall from the example depicted in 
Figure 1 that 2 P is a lower bound for polynomials of degree p> 2. 





maximum latency 


average latency 


unfairness 


min-max flow 


1 


a{C) 


? 


system optimum 


7(£) 


1 


7(£) 


user equilibrium 


a(£) 


a{C) 


1 
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Abstract. The problem of finding the value and optimal strategies of 
players in cyclic games is studied. A polynomial time algorithm for solv- 
ing cyclic games is proposed. 



1 Introduction and Problem Formulation 



Let G = (V, E) be a finite directed graph in which every vertex u £ V has 
at least one leaving edge e = (u,v) € E. On edge set E is given a function 
c: E — » R which assign a cost c(e) to each edge e £ E. In addition the vertex set 
V is divided into two disjoint subsets V A and V B (V = V A U V B , V A fl V B = 0) 
which we will regard as a positions sets of two players. 

On G we consider the following two-person game from [1,2]. The game start 
at position Vq £ V. If vq £ V A then the move is done by first player, otherwise 
it is done by second one. The move means the passage from position Vq to the 
neighbour position Vi through the edge ei = (vo,Vi) £ E. After that if Vi £ V A 
then the move is done by first player, otherwise it is done by second one and 

so on indefinitely. The first player has the aim to maximize lim inf — c(ej) 



while the second player has the aim to minimize lim sup — c(ei). 

t—> OO t ' ^ 



i= 1 



In [1,2] is proved that for this game there exists a value p(v o) such that the 

first player has a strategy of moves that insures lim inf — N c(ej) > p(vo) and 

t— >oo t 

i—1 

the second player has a strategy of moves that insures lim sup - } c(ej) < 

t — >oo t 

i = 1 

p{v o). Furthermore in [1,2] is shown that the players can achieve the value p(v o) 
applying the strategies of moves which do not depend on t. This means that the 
considered game can be formulated in the terms of stationary strategies. Such 
statement of the game in [2] is named cyclic game. 

The strategies of players in cyclic game are defined as a maps 



: u 



■ v £ V G (u) for u £ V A ; s B : u — > v £ V G (u) for u £ V B 



D. Bienstock and G. Nemhauser (Eds.): IPCO 2004, LNCS 3064, pp. 74—85, 2004. 
(c) Springer- Verlag Berlin Heidelberg 2004 




Polynomial Time Algorithm 



75 



where V G (u) represents the set of extremities of edges e = (u, v ) £ E, i.e. V G (u) = 
{»£ V | e = (u, v) £ E}. Since G is a finite graph then the sets of strategies of 
players 

S A = {s A : u — > v £ V G (u) for u £ V A }; S B = {s B : u — > v £ V G (u) for u £ V B } 

are finite sets. The payoff function F Vo : S A x S B — > R in cyclic game is defined 
as follows. 

Let s A £ S A and s B £ S B be a fixed strategies of players. Denote by G s = 
(V, E s ) the subgraph of G generated by edges of form (u, s A (u)) for u £ V A and 
(u,s B (u)) for u £ V B . Then G s contains a unique directed cycle C s which can 
be reached from vo through the edges e £ E s . The value F Vq (s a , s b ) we consider 
equal to mean edges cost of cycle C s , i.e. 



Fvq (Sa? $ B ) 



1 

n(C s ) 



c ( e )> 

e£E(Cs) 



where E(C S ) represents the set of edges of cycle C s and n{C s ) is a number of 
the edges of C s . So, the cyclic games is determined uniquely by the network 
(G, V A .V B . c) and starting position vq. In [1,2] is proved that there exist the 
strategies s* A £ S A and s% £ S B such that 



p(v) = F v (s* a , s%) = max min F v (s A ,s B )= min max F v (s A ,s B ), Mv£V. 
s a &S a s b gS b s b GS b s a gS a 

So, the optimal strategies s A ,s* B of players in cyclic games do not depend on 
starting position v although for different positions u. v £ V the values p(u) 
and p(v) may be different. It means that the positions set V can be divided 
into several classes V = V 1 U V 2 U ... U V k according to values of positions 
p 1 ,p 2 1 ■ ■ ■ ,p k , i.e. u,v £ V 1 if and only if p l = p{u) = p(v). In the case k = 1 
the network (G. V A ,V B . c) is named the ergodic network [2]. In [3] is shown 
that every cyclic game with arbitrary network (G. V A .V B . c) and given starting 
position Vn can be reduced to an auxiliary cyclic game on auxiliary ergodic 
network (G',V A ,V B ,c'). 

It is well-known [2,4] that the decision problem associated to cyclic game 
is in NPnco-NP but it is not yet known for this problem to be in P. Some 
exponential and pseudo-exponential algorithms for finding the value and the 
optimal strategies of players in cyclic are proposed in [1,2,4]. 

A similar problems on acyclic networks have been considered in [5,6]. A 
new classes of acyclic games have been studied in [5,6] and polynomial time 
algorithms for its solving have been elaborated. Moreover in [5,6] is made an 
attempt to reduce the cyclic game to acyclic one. Unfortunately the reduction 
from [5,6] gives the first player some control over the length of the formed cycles, 
so that it is not conservative in the general case. 

In this paper we propose a polynomial time algorithm for finding the value 
and optimal strategies of players in cyclic games. The proposed algorithm is 
based on a new reduction from cyclic games to acyclic one. 
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2 Some Preliminary Results 

First of all we need to remind some preliminary results from [2]. 

Let (G, V A , V B ,c) be a network with the properties described in section 1. We 
denote 

{ max {c(m,u)} for u £ V A , 
vev G (u) 

min {c(u, u)} for u£V B , 
i.GVgG) 

VEXT(c,u ) = {u £ V G (u) | c(u,v) = ext(c,u)}. 

We shall use the potential transformation d(u , v) = c(u, v) + s(v) — s(u) for costs 
on edges e = (u, v) £ E, where e: V — * R is an arbitrary function on vertex set 
V. In [2] is noted that the potential transformation do not change the value and 
the optimal strategies of players in cyclic games. 

Theorem 1. Let (G,V A ,V B ,c) be an arbitrary network with the properties de- 
scribed in section 1. Then there exists the value p(v), v £ V and the function 
e: V — > R which determine a potential transformation c'(u,v) = c(u,v ) +e(f) — 
e(u) for costs on edges e = (i i,v) £ E, such that the following properties holds 

a) p(u) = ext{c ' , u) for v£V, 

b) p\u) = p(v) for u£V a UV b and v £ VEXTfd ,u), 

c) p(u) > p(v) for u £ V A and v £ V G (u), 

d) p(u) < p(v) for u £ V B and v £ V G (u), 

e) max|c'(e)| < 2|V| max |c(e)| . 

e£E e£E 

The values p(v),v £ V on network (G,V A ,V B ,c) are determined unequally and 
the optimal strategies of players can be found in the following way: fix the arbi- 
trary strategies s* A : V A — > V and s* B : V B — » V such that s* A (u) £ VEXT{d ,u) 
for u £ V A and s* B (u) £ VEXT(c',u ) for u £ V B . 

The proof of Theorem 1 is given in [2] . 

Further we shall use the theorem 1 in the case of the ergodic network 
(G, Vi,V 2 ,c), i.e. we shall use the following corollary. 

Corollary 2. Let (G, V A .V B . c) be an ergodic network. Then there exists the 
value p and the function e: V — ■> R which determine a potential transformation 
c'(u,v) = c(u,v) + e(v) — e(u) for costs of edges e = (u,v) £ E such that p = 
ext{c' , u) for u £ V. The optimal strategies of players can be found as follows: fix 
arbitrary strategies s A : V A —> V and s* B : V B — > V such that s* A (u) £ VEXT(c', u) 
for u £ V A and s* B (u ) £ VEXT(cf ,u) for u £ V B . 

3 The Reduction of Cyclic Games to Ergodic Ones 

Let us consider an arbitrary network (G. V A ,V B , c) with given starting position 
Vo £ V which determine a cyclic game. In [3] is shown that this game can be 
reduced to a cyclic game on auxiliary ergodic network (G\ W A , W B , c), G' = 
(' W , F) in which is preserving the value p(v o) and vq £ W = V U X U Y. 
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The graph G' = ( W , F) is obtained from G if each edge e = (u, v ) is changed 
by a triple of edges e 1 = (u,x),e 2 = ( x,y),e 3 = {y,v) with the costs c^e 1 ) = 
c(e 2 ) = c(e 3 ) = c(e). Here x £ G Y and u,v £ V;W = V U X U Y. In 
addition in G' each vertex x is connected with vo by edge ( x , Vq) with the cost 
c(x,v o) = M ( M is a great value) and each edge (y,v o) is connected with vo by 
edge (y, vo) with the cost c = (y, vo) = —M. In (G' , W A , W B ,c) the sets W A and 
W B are defined as follows: W A = V A UY; W B = V B U X. 

If easy to observe that this reduction can be done in linear time. 

4 Acyclic Games on Networks and Polynomial 
Algorithms for Its Solving 

We consider the following auxiliary games from [5,6]. 

4.1 Acyclic c-Game on Acyclic Networks 

Let G = (V,E) be a finite directed graph without directed cycles. Assume that 
in G there exists a vertex v f which is attainable from each vertex v £ V, i.e. v f 
is a sink in G. In addition we consider that on edge set E is given a function 
c: E — > R and the vertex set V is divided into two disjoint subsets V A and 
V B (V = Va U Vb, V a l~l V B = 0) which we regard as positions sets of two players. 

On G we consider the two-person game from [5]. The game start at given 
position vo G P\ {u r }. If vq G V a then the move is done by first player otherwise 
it is done by second one. Like in cyclic game here the move means the passage 
from position vo to the neighbor position v\ through the edge ei = (uo,Ui). 
Again if v A G V A then the move is done by first player, otherwise it is done 
by second one an so on while the final position v f is not reached. As soon as 
the final position v f is reached the game is over. The final step of the game we 
denote by t. In this game the first player has the aim to maximize the integral 

t ' t 

cost ^^c(ei) and the second player has the aim to minimize ^^c(e,). 

i—l i=l 

In [5] is shown that for such game there exists a value p(v o) such that the 

t 

first player has a strategy of moves that insures ^ c(e,) > p(v o) and the second 

»= 1 
t 

player has a strategy of moves that insures ^c(ei) < p(v o). Moreover the 

i — 1 

players can active this values using the stationary strategies of moves. This 
game in [5] is named acyclic c-game. The strategies of players in acyclic c-game 
can be defined as a maps 

s A : u — » v G V G {u) for u G V A \ {uo}; s B : u — + v G V G (u) for u G V B \ {r’o}- 
Since G is a finite graph then the set of strategies of players 

S A = u > v £ V G (u) for u £ V A \ {u 0 }}; 
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S B = {s B : «-*»£ V G (u) for u £ V B \ {v 0 }} 

are finite sets. The payoff function F Vo : S A x S B —> R in acyclic c-game is defined 
as follows. 

Let s A € S A and s B € S B be a fixed strategies of players. Denote by T s = 
(V. E s ) the tree with root vertex v f in G generated by edges of form (u, s A (u)) 
for u € V A \ {z>o} and (u,s B (u)) for u £ V B \ {fo}- Then in T s there exists a 
unique directed path P Ta (v o,v f ) from vq to v f . We put 

F Vo {s a ,s b )= ^2 c ( e )’ 

eGE^Pr,, (vo 

where E(P Ts ( vo , v f )) represents the set of edges of the directed path P Ts (vo, v f ). 

In [5,6] is proved that for acyclic c-game on network ( G , V\, V 2 , c) there exist 
the strategies of players s*,s* such that 

p(v) = F v {s* a , s* b ) = max min F v (s A ,s B ) = min max F v (s A ,s B ) (1) 
s A £S a s b gS b s b gSb s A eS A 

and s*,s* do not depend on starting position v £ V, i.e. (1) holds for every 
v £ V . Note that here for different positions u and v the values p(u) and p(v) 
may be different. 

The equality (1) is evident in the case when ext(c, u) = 0, V u £ V. In this 
case p(u) = 0, V u £ V and the optimal strategies of players can be obtained by 
fixing the maps s*: V A \ {17} — » V and s* : V B \ {r^} — » V such that s* (u) £ 
VEXT[c,u ) for u £ V A \ {r> / } and s%(u ) £ VEXT{c,u ) for u £ V B \ {c/}. 

If the network (G,V A ,V B ,c) has the property ext(c,u ) = 0, Vm £ V \ {r> / } 
then it is named the network in canonic form [5]. So, for the acyclic c-game on 
network in canonic form the equality (1) holds and p(v) = 0, Vw £ V. 

In general case the equality (1) can be proved using the properties of the 
potential transformation c'(u,v) = c(u,v ) + e(v) — e(u) on edges e = (u,v) of 
the network, where s: V — > R is an arbitrary real function on V . The fact is 
that such transformation of the costs on edges of the network in acyclic c-game 
do not change the optimal strategies of players, although the values p(v) of the 
positions v £ V are changed by p(v) + s(vf) — s(v). It means that for arbitrary 
function e: V — » R the optimal strategies of the players in acyclic c-games on 
networks (G, V A , V B ,c ) and on (G, V A , V B . c') are the same. Using such property 
in [5] is proved the following theorem. 

Theorem 3. For arbitrary acyclic network (G,V A ,V B ,c) there exists a function 
e: R — > R which determine a potential transformation c'(u, v ) = c(u, v) + e(v) — 
e(u) on edges e = (u, v) £ E such that (G, V A ,V B ,c') has the canonic form. 
The value e(v),v £ V, which determine e: V — > R, can be found by using the 
following recursive formula 

( e{v f ) = 0, 

J ( max {c(u,v) +£(«)}, for u £ V A \ {v f }, 

\ £ ( u ) = l v&V G (u) 1 ) 

^ j min {c(u, v) + e(v)}, for u £V B \{v f }. 

\ I i 'u6V g (m) 
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According to this theorem the optimal strategies of players in acyclic c-game 
can be found using the following algorithm. 

Algorithm 1. 

1. Find the values e(v),v £ V, according to recursive formula (2) and the 
corresponding potential transformation c'(u, v ) = c(u, v) + e(v) — e(u) for edges 
(u, v ) £ E. 

2. Fix the arbitrary maps s* A : V A \ {v f } — > V and s%: V B \ {t) f } — * V for 
which s A (u) € VEXT(c',u) for u € V A \ {v f } and s* B (u) £ VEXT(c',u) for 
u£V b \ E}. 

Remark. The values e(u),u £ V, represent the values of the acyclic c-game on 
(G,V a ,V b ,c) with starting position u, i.e. e(u) = p(u), V« G V. Algorithm 1 
need 0(n 2 ) elementary operations because such number of operations need the 
tabulation of values e(u),u G V, by using formula (2). 

4.2 Acyclic c-Game on Arbitrary Networks 

Let us consider the game from section 4.2 when G is an arbitrary graph, i.e. G 
may contains directed cycles. In this case the subgraph T s = (V,F S ), generated 
by given strategies of players 

s A : u — > v G V G (u) for u € V A \ {«/}; s B : u —> v G V G (u) for u G V B \ {r’/} ; 

may contains directed cycles. So, for given s A and s B either exists a unique 
directed path P Ts (vo, v f ) from vq to v f in T s or such a path does not exist in T s . 
In the second case there exists a unique directed cycle C s which can be reached 
from Vo through the edges of T s . If in T s there exist a unique directed path 
Pt s (v o,v f ) from vq to v f then F Vo (s A , s B ) we define as 

F Vo (s a ,s b ) = ^2 c ( e )- 

e€E^PT s (vo < v f)) 

If in T s there are no directed paths from vq to v f then F Vo (s A , s B ) we define in 
the following way. 

Let C s be the cycle which can be reached from Vo and P Ts (vo, uq) represent a 
unique directed path in T s which connect vo with the cycle C s (P T (vo,Uo) and 
C s have no common edges) . Then we put 

+00 if E c(e) > 0, 

eeE{C s ) 

E c ( e ) if E c ( e ) = °> 

eSi5(P^ s (vo,u 0 )) eeE(C s ) 

—oo if ^2 c ( e ) < 0- 

eeE(Cs) 
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We consider the problem of finding the strategies s* , s* for which 

Fy 0 (s* A ,s* B ) = max min F V0 (s A ,s B ). 
sa€Sa sb€Sb 

This problem [7] is named the max-min path problem in c-game. 

It is easy to observe that for this problem may be the case when the players 
couldn’t reach the final position v f . We shall formulate the condition when for 
the player in this game there exists a settle value p(v) for each v £ V. 

Theorem 4. Let (G,V A ,V B ,c) be a network with given final position, where G 
is on arbitrary directed graph. Moreover let us consider that c(e) ^ 0 for 

e<EE{C s ) 

every cycle C s from G. Then for c-game on (G,V A ,V B ,c) the conditions 

p(v) = F v (s* A , s* B ) = max min F v (s A ,s B ) = min max F v (s a ,s b ),\/vGV 
s a GS a sb£Sb sb&Sb sa&S a 

holds if and only if there exists a function e: V — » R which determine a po- 
tential transformation c'{u,v ) = c(u,v) + e(v) — e(u) on edges (u,v) £ E such 
that ext(c!,u) = 0,\/v € V. The optimal strategies of players in c-game on 
( G , V A , V B ,c) can be found by fixing the arbitrary maps V A \ { 17 } — > V and 
s* B : V B \ {v f } — y V such that s A (u) £ VEXT(c', u) for u £ V A and s B (v ) £ 
VEXT(c',u) for u £ V B . 

Proof. =>• Let us consider that the condition 

p(v) = F v (s*.,s* R ) = max min F v (s A ,s B ) = min max F v (s A ,s B ) 
s A &S a sb&Sb sb&Sbs a gS a 

holds for every v £ V and p(v) < 00 , Vu £ V. It is easy to observe that if we 
put e(v) = p(v) for v £ V then the function e: V — * R determine a potential 
transformation c'(u, v ) = c(u, v)+£(v) — s(u) on edges ( u , v) £ E such that holds 
ext(c', u) = 0, Vu £ V. 

4= Let us consider that there exists a potential transformation c' (u, v) = 
c(u,v)+s(v)—e(u) of costs of edges (u,v) £ E such that ext(c',u) = 0, V« £ V. 
Then for a c'-game on network (G, V A ,V B , c') the following condition 

max min F'(s A , s B ) = min max F'(s A , s B ), Wv £ V 
s a GS a sb&Sb s b &Sbs a £Sa 

holds. Since the potential transformation do not change the optimal strategies 
of players then 

max min F v (s A ,s B )= min max F v (s A , s B ), Vu £ V. 

sa&S a sb&Sb sb&Sb sa&Sa 



Corollary 5. The difference e(v) — e(vf),v € V, in Theorem 4 determine the 
cost of max-min path from v to v f , i.e. e(v) = p(v),\/v £ V. 
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Let us consider the network (G, V A , V B , c ) for which there exist min-max 
paths from every v £ V to given Vg £ V. Then the following algorithm finds all 
values s(v), v € V, of min-max paths from v to Vg. 

Algorithm 2. 

1. We construct an auxiliary acyclic graph H = (IT, E), where 

W = {wg} u w 1 u w 2 u . . . u w n , w i n W j = 0, t 7^ j; 

w* = {w z 0 ,w\,...,w l n _ J, i = l,n; 

E = E° U E 1 U E 2 U . . . U E n ~ l - 

E z = {(«4 +1 ,u;|)| (vk,vi) e E}, i = l,n - 1; 

E ° = {KU w o)l (.v k ,v 0 ) € E, i = 1, n}. 

The vertex set IT of H is obtained from V if it is doubled n times and then a 
sink vertex u>g is added. The edge subset E l C E in H connect the vertices of 
the set 1T* +1 and the vertices of the set IT* in the following way: if in G there 
exists an edge (vk,vi) € E then in H we add the edge The edge 

subset E° C E in H connect the vertices w\ £ W 1 U IT 2 U . . . U W n with the 
sink vertex wj] , i.e. if there exists an edge (vk,i’o) €E E then in H we add the 
edges (w l k ,Wg) £ E°, i = l,n. 

After that we define the acyclic network (Ho, W A , W B , Co), H 0 = (Wo, So) 
where Ho is obtained from H by deleting the vertices w l k £ W from which 
the vertex u;}] couldn’t be attainable. The sets W A ,W B and the cost function 
Co: E 0 — > R are defined as follows: 

W A = {wi £ Wo | v k £ V A }, W B = {wi £ Wo | v k £ Tb}; 

co(w l k hl ,w k ) = c(v k ,vi) if (v k ,vi) £ E and (w^ 1 , w z k ) £ E 0 ; 
c 0 (w l k ,Wg) = c(v k ,v 0 ) if (v k ,v 0 )€E and (w k ,Wg) £ E 0 . 

2. On auxiliary acyclic network (Ho, W A , W B , Co) with sink vertex icq we 
consider the acyclic c-game. According to recursive formula (2) we find s(w l k ) 
for every w\ £ H 0 , i.e. e(w k ) = p(w l k ), V w k £ H 0 . 

3. Fix U = {t'o}; i = 0. 

4. i = i+l. 

5. Find the edge set E(U) = {(v k ,vi) £ E\vi £ U,v k £ V\U} and the 
vertex set V(U) = {v k £V\e= (v k ,vi) £ E(U)}. 

6. Fix e(v k ) = e(w k ) for v k £ V(U). 

7. Select the vertex set V°(U) from V(U) where v k £ V°(U) if the following 
conditions are satisfied: 

a) c(v k ,v q ) + e(v q ) - e(v k ) < 0, \/ v q £ U U V(U), (v k ,v q ) £ E if v k £ 

v A nv(uy, 

b) c(v k ,v q ) + e(v q ) - e(v k ) > 0, M v q £ U U V(U), (v k ,v q ) £ E if v k £ 

v B nv(uy, 

c) there exist an edge (v k , Vi) € E such that c(v k , Vi)-\-£(vi)—£(v k ) = p(v o). 
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8. If V°(U) = 0 go to step 10. 

9. Change U by U U V°(U). 

10. If U ^ V then go to step 4, otherwise go to next step 11. 

11. End. 

The correctness of this algorithm follow from the algorithm from [8] . 



4.3 Acyclic i-Game on Networks 

Let (G, V A ,V B , c) be the network with the properties described in section 1. So, 
G = ( V , E ) represent a directed graph without directed cycles and G contains 
a sink vertex Vf £ V. On E is defined a function c: E — > R and on V is given 
a partition V = V A U V B (V A fl V B = 0) where V A and V B are considered as a 
positions sets of two player. 

We study the following acyclic l- game from [5] . Again the strategies of players 
we define as a maps 

s A : u — » v € V G (u) for u £ V A \ {u/}; s B : u —> v £ V G (u) for u £ V B \ {u/}- 

The payoff function F Vo : S A x S B — > R in this game we define as follows. 

Let s A £ S A and s B £ S B be a fixed strategies of players. Then the tree 
T s = (V,E S ) contains a unique directed path P Tb (v o,v) from vq to v f . We put 

F ”"(»—) = „ ( P J m , v , n S 

eeE(P Ts (v 0 ,v f )) 

where E ( P Ts (uo, v f )) is the set of edges of the path P Ts (vq, v f ) and n ( P Ta (vo, v f )) 
is the number of its edges. 

In [5,6] is proved that the function F Vo (s A , s B ) satisfy the following condition: 

¥(v o) = Fv 0 (s2, s*) = max SA6 s A min SBe s B F Vo (s A ,s A ) = 

w) 

= min SB6 5 B ma x saGSa F Vo (s a ,s b ). 

In this game the optimal strategies sj[ and s B of players depend on starting 
position vq- But the players can achieve the settle value p(v o) starting from vq 
by using the successive choosing of optimal strategies of moves from one position 
to another while the final position v f is not reached. 

The equality (3) can be proved by using the Theorem 1 (corollary 2). It is 
easy to observe that if for given acyclic l- game we identify the positions Vq and 
v f we obtain the cyclic game with an auxiliary ergodic network (G', V A , V' B1 d) 
and given starting position vo- In this ergodic network the graph G' is obtained 
from G when the vertices vq and v f are identified. In the graph G' all directed 
cycles pass through the vertex vq- Applying the corollary 1 of Theorem 1 to the 
obtained cyclic games on ergodic network (G', V' A -V' B . c!) with starting position 
Vq, we obtain the proof of the following theorem [5]. 
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Theorem 6 . Let be given an acyclic l -game determined by the network 
(G, V A ,V B , c) with starting position vo- Then there exist a value p(v o) and a func- 
tion e: V — » R which determine a potential transformation c'(u,v) = c(u,v) + 
e(v) — e(u) of costs on edges e = (u,v) € E such that the following conditions 
holds 

a) p(i’o) = ext(c',u), Vu € V \ {v f }; 

b) e(v 0 ) = e(v,) 

the optimal strategies of players in acyclic l-game can be found as follows: fix 
the arbitrary maps V A \ {v/} — > V and s* B : V B \ {v f } — > V such that s A (u) £ 
VEXT(c’ ,u) for u£V \ { 17 } and s B (u) £ VEXTfd ,u). 

This theorem follow from corollary 2 of Theorem 1 because the acyclic /- 
game on network Let (G, V A ,V B . c:) with starting position vq can be regard as 
a cyclic game on network Let (G', V A , V B , d) with starting position vq. This 
network is obtained from the network (G, V A ,V B , c) when the vertices vq and v f 
are identified. 

Taking in consideration remark to Algorithm 1 and corollary 5 of Theorem 4 
we obtain the following result. 

Corollary 7. The difference e(v) — e(vo),v £ V, in Theorem 6 represent the 
costs of max-min path from v to v f in acyclic c-game on network (G, V A , V B ,c) 
where c(u, v ) = c(u, v) — pf vq), V ( u , v) £ E. 

So, if p(v 0 ) is known then the optimal strategies and s B of players in 
acyclic /-games can be found in time 0 (n 2 ) by reducing acyclic /-game to acyclic 
c-game on (G, V A ,V B ,c). 

On the basis of this property in [5] is elaborated polynomial algorithm for 
solving /-games. In [6] is made an atempt to elaborate strongly polynomial time 
algorithm for /-games. Unfortunately the algorithm from [6] is not valid. 

5 Polynomial Time Algorithm for Solving Cyclic Games 

On the basis of the obtained results we can propose polynomial time algorithm 
for solving cyclic games. 

We consider an acyclic game on ergodic network (G,V A ,V B ,c) with given 
starting position Vq . The graph G = (V,E) is considered to be strongly con- 
nected and V = {vo,Vi,V 2 , ■ ■ ■ ,v n -i}. Assume that vq belong to the cycle G s * 
determined by the optimal strategies of players s* and s* . 

We construct an auxiliary acyclic graph H r = ( W r ,E r ), where 

W r = {wg} U W 1 U W 2 U . . . U W r , W i nW j = 0, i^ j; 

i=l~r; 

E = E° U E 1 U E 2 U . . . U E r ~ 1 -, 

Ei = {(wl +1 ,wi)\(v k ,vi) £ E}, i = l,r - 1; 

E ° = (K^o)l ( v k,v 0 ) e E, i = l~F}. 
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The vertex set W r of H r is obtained from V if it is doubled r times and then a 
sink vertex Wq is added. The edge subset E l C E in H r connect the vertices of 
the set W l+1 and the vertices of the set W l in the following way: if in G there 
exists an edge (vk,vi) £ E then in H r we add the edge ( w ,w\). The edge 
subset E° C E in H r connect the vertices w\ £ W 1 U W 2 U . . . U W r with the 
sink vertex u>q, i.e. if there exists an edge (vk,vo) £ E then in H r we add the 
edges (wl,Wg) £ E°, i = 1 , r. 

After that we define the acyclic network (H' r , W A , W B , d r ), H' r = ( W r ,E r ) 
where H' r is obtained from H r by deleting the vertices w k £ W r from which 
the vertex Wq couldn’t be attainable. The sets W A ,W B and the cost function 
c' r : E r — > R are defined as follows: 

W A = {wi £ Wo | v k £ V A }, W B = {wi £ Wo I V k £ V B y, 
d r {w l ^, w l k) = c{vk,vi) if (vk,v t ) £ E and w l k ) £ E r \ 

c' r (wl,Wo) = c(v k ,v 0 ) if (vk,v 0 ) £ E and (w k ,Wo) £ E r . 

Now we consider the acyclic c-game on acyclic network (H' r , W A . W B ,c' r ) with 
sink vertex Wq and starting position Wq. 

Lemma 8. Let p = p(v o) is a value of the ergodic cyclic game on G and the 
number of edges of the max-min cycle in G is equal r. Moreover, let p,.(Wq) 
be the value of the l-game on (H' r ,W A ,W B ,c' r ) with starting position Wq. Then 

P(V o) =Pr(w o). 

Proof. It is evident that there exist a bijective mapping between the set of 
cycles with exactly r edges (which contains the vertex uq) in G and the set of 
directed paths with r edges from Wq to w § in H' r . Therefore p(v o) = p t {wq). 

On the basis of this lemma we can propose the following algorithm for finding 
the optimal strategies of players in cyclic games. 

Algorithm 3. 

We construct the acyclic networks ( H ' r , W A , W B ,c' r ), r = 2, 3, . . . , n and for 
each of them solve Z-game. In a such way we find the values P 2 {wq),ps(wq), . . . , 
p n (w q) for these Z-games. Then we consequtively fix p = P 2 (wq) , P 3 (wq ) , . . . , 
Pn{w o) and each time solve the c-game on network (G, V A , V B , d), where c' = 
c — p. Fixing each time the values d{vk) = jf{vk) for vu £ V we check if the 
following condition 

ext(c',Vk) = 0. \/ Vk £ V 

is satisfied, where c'{vk,vi) = d(vk,vi) + e{vi) — e(vk)- We find r for which 
this condition holds and fix the respective maps and such that s* A {vk) £ 
VEXT(d, Vk) for u £ V A and s* B (vk) £ VEXT(c',Vk ) for u £ V B . So, s* A and s* B 
represent the optimal strategies of player in cyclic games on G. 

Remark. Algorithm 3 find the value p(v o) and optimal strategies of players in 
time 0(n 2 \E\ 2 ). 
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Abstract. We propose a general methodology based on robust opti- 
mization to address the problem of optimally controlling a supply chain 
subject to stochastic demand in discrete time. The attractive features 
of the proposed approach are: (a) It incorporates a wide variety of phe- 
nomena, including demands that are not identically distributed over time 
and capacity on the echelons and links; (b) it uses very little information 
on the demand distributions; (c) it leads to qualitatively similar optimal 
policies (basestock policies) as in dynamic programming; (d) it is numer- 
ically tractable for large scale supply chain problems even in networks, 
where dynamic programming methods face serious dimensionality prob- 
lems; (e) in preliminary computational experiments, it often outperforms 
dynamic programming based solutions for a wide range of parameters. 



1 Introduction 

Optimal supply chain management has been extensively studied in the past us- 
ing dynamic programming, which leads to insightful policies for simpler systems 
(basestock policies for series systems; Clark and Scarf [7]). Unfortunately, dy- 
namic programming assumes complete knowledge of the probability distributions 
and suffers from the curse of dimensionality. As a result, preference for imple- 
mentation purposes is given to more intuitive policies that are much easier to 
compute, but also suboptimal (see Zipkin [12]). 

Hence, the need arises to develop a new optimization approach that incor- 
porates the stochastic character of the demand in the supply chain without 
making any assumptions on its distribution, is applicable to a wide range of 
network topologies, is easy to understand intuitively, and combines computa- 
tional tractability with the structural properties of the optimal policy. The goal 
of this paper is to present such an approach, based on robust linear and mixed 
integer optimization that has witnessed increased research activity (Soyster [11], 
Ben-Tal and Nemirovski ([1,2,3]) and El-Ghaoui et. al. ([8,9], Bertsimas and Sim 
[5,6]). We utilize the approach in [5,6], which leads to linear robust counterparts 
while controlling the level of conservativeness of the solution. 

The contributions of this paper are as follows: (a) We develop an approach 
that incorporates demand uncertainty in a deterministic manner, remains nu- 
merically tractable as the dimension of the problem increases and leads to high- 
quality solutions without assuming a specific demand distribution, (b) The ro- 
bust problem is of the same class as the nominal problem, that is, a linear 
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programming problem if there are no fixed costs or a mixed integer program- 
ming problem if fixed costs are present, independently of the topology of the 
network, (c) The optimal robust policy is qualitatively similar to the optimal 
policy obtained by dynamic programming when known. In particular, it remains 
basestock when the optimal stochastic policy is basestock, as well as in some 
other cases where the optimal stochastic policy is not known, (d) We derive 
closed-form expressions of key parameters defining the optimal policy. These ex- 
pressions provide a deeper insight into the way uncertainty affects the optimal 
policy in supply chain problems. 



2 The Robust Optimization Approach 

We rely extensively on the robust optimization tools developed by Bertsimas and 
Sim in [5] for linear programming problems. We consider the following problem 
subject to data uncertainty: 

min c'x : Ax < b, 1 < x < u, 

where we assume WLOG that only the matrix A is subject to data uncertainty. 
Let A = {A G K mxn | a,ij G [<% - % , a ^ + a ] \/i,j, ^ \ ai i ^ a dl < r}. 

_ (iJ)eJ Uij 

r is a parameter that controls the degree of conservatism. The robust problem 
is then formulated as: 

min c'x 

s.t. Ax < b, VA G A (1) 

1 < x < u. 



Theorem 1 (Bertsimas and Sim [5]). The uncertain linear programming 
problem has the following robust, linear counterpart. : 
min c'x 

s.t. ^ hijXj + qiT + ^ nj < bi, Vi 

qi + nj > a,ijyj, V(i, j) G J (2) 

-y < x < y 

1 < X < u 

q > 0, r > 0, y > 0. 

The robust counterpart is therefore of the same class as the nominal prob- 
lem, that is, a linear programming problem. This is a highly attractive feature 
of this approach, since linear programming problems are readily solved by stan- 
dard optimization packages. Moreover, if in the original problem (1), some of 
the variables were constrained to be integers, then the robust counterpart (2) 
would retain the same properties, i.e. , the robust counterpart of a mixed integer 
programming problem is itself another mixed integer programming problem. 
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3 The Single Station Case 



3.1 The Uncapacitated Model 

In this section we apply the robust optimization framework to the problem of 
ordering, at a single installation, a single type of item subject to stochastic 
demand over a finite discrete horizon of T periods, so as to minimize a given 
cost function. We define, for k = 0, . . . ,T: 

Xk '■ the stock available at the beginning of the kth period, 

Uk : the stock ordered at the beginning of the kth period, 

Wk : the demand during the fcth period. 

The stock ordered at the beginning of the fcth period is delivered before the 
beginning of the (k + l)st period, that is, all orders have a constant leadtime 
equal to 0. Excess demand is backlogged. Therefore, the evolution of the stock 
over time is described by the following linear equation: 

x k +i = x k + u k - w k , k = 0, . . . ,T — 1, (3) 



leading to the closed- form expression: 



k 

X k + 1 = X 0 + - Wi), k = 0,#..,T-l. (4) 

i= 0 

Neither the stock available nor the quantity ordered at each period are subject 
to upper bounds. Section 3.2 deals with the capacitated case. 

The demands w k are random variables. In order to apply the approach out- 
lined in Sect. 2, we model w k for each k as an uncertain parameter that takes 
values in [Wfc — Wk,w k + u>k\- We define the scaled deviation of w k from its nom- 
inal value to be z k = (w k — Wk)/wk, which takes values in [—1,1], We impose 
budgets of uncertainty at each time period k for the scaled deviations up to 
time k. Hence, we now have the constraint y~^ fc _ 0 \z- t \ < for all time periods 
k = 0 , ,T — 1. These budgets of uncertainty rule out large deviations in the 
cumulative demand, and as a result the robust methodology can be understood 
as a “reasonable worst-case” approach. The main assumption we make on the r k 
is that they are increasing in k, i.e. , we feel that uncertainty increases with the 
number of time periods considered. We also constrain the to be increasing 
by at most 1 at each time period, i.e., the increase of the budgets of uncertainty 
should not exceed the number of new parameters added at each time period. 

Finally, we specify the cost function. The cost incurred at period k consists 
of two parts: a purchasing cost C(u k ) and a holding/shortage cost resulting from 
this order R(xk+ 1 )- Here, we consider a purchasing cost of the form: 



f I\ + c ■ u, if u > 0, 
\ 0, if u = 0, 



( 5 ) 



with c > 0 the unit variable cost and K > 0 the fixed cost. If I\ > 0, a 
fixed positive cost is incurred whenever an order is made. The holding/shortage 
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cost represents the cost associated with having either excess inventory (positive 
stock) or unfilled demand (negative stock). We consider a convex, piecewise 
linear holding/shortage cost: 

R(x) = max(hx, — px), (6) 



where h and p are nonnegative. The holding/shortage cost for period k, yk, is 
computed at the end of the period, after the shipment Uk has been received 
and the demand Wk has been realized. We assume p > c, so that ordering stock 
remains a possibility up to the last period. 

Using the piecewise linearity and convexity of the holding/shortage cost func- 
tion, and modelling the fixed ordering cost with binary variables, the inventory 
problem we consider can be written as a mixed integer programming problem: 



T-l 



min E (cu k + Kv k + Vk) 

k = 0 

s.t. y k >h(x o + y - ic)^ k = 0, . . .,T - 1 



*= o 



( 7 ) 



Vk > -p ( x 0 + - m) ] k = 0, , . . , T — 1 



\ i—0 / 

0 < u k < Mv k , Vk € {0, 1} k = 0, . . . , T - 1, 

where Wi = vJi + Wi- Zi such that z £ V = {\zi\ < 1 Vi > 0, ^^=0 \ Zi \ — Tfe Vk > 
0}. Applying Theorem 1, we obtain: 

Theorem 2. The robust formulation for the single- station inventory problem 
(7) is: 

T-l 

min E (' cu k + Kv k + Vk) 

k—0 

/ k k \ 

s.t. Vk>h\x o + y^(u* - Wi) + q k r k + y: r ik 



i — 0 



i= 0 



fc k \ 

yk > pi -X 0 - y ^(u* - ttj») + Qk r k + r^ 

\ i—0 i = 0 y 

+ r ik > Wi 

qk > 0, r ik > 0 
0 < u k < Mv k , v k G {0, 1}, 



(8) 



where M is a large positive number. 



The variables qk and r ik quantify the sensitivity of the cost to infinitesimal 
changes in the key parameters of the robust approach, specifically the level of 
conservativeness and the bounds of the uncertain variables, qk Tk + r *fc 

represents the extra inventory (or lack thereof) that we want to take into account 
in controlling the system from a worst-case perspective. 
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The robust problem is a linear programming problem if there is no fixed cost 
( K = 0) and a mixed integer programming problem if fixed costs are present 
( K > 0). In both cases, this robust model can readily be solved numerically 
through standard optimization tools, which is of course very appealing. It is 
also desirable to have some theoretical understanding of the optimal policy, in 
particular with respect to the optimal nominal policy and, if known, the optimal 
stochastic policy. We address these questions next. 

Definition 1 ((S,S) and (s,S) policies). The optimal policy of a discrete- 
horizon inventory problem is said to be (s,S), or basestock, if there exists a 
threshold sequence ( Sk,Sk ) such that, at each time period k, it is optimal to 
order Sk — Xk if Xk < Sk and 0 otherwise, with Sk < Sk ■ If there is no fixed 
ordering cost (I\ =0), Sk = Sk- 
in order to analyze the optimal robust policy, we need the following lemma: 

Lemma 1. (a) The optimal policy in the stochastic case, where the cost to min- 
imize is the expected value of the cost function over the random variables Wk, is 
( s,S ). As a result, the optimal policy for the nominal problem is also (s,S). 

(b) For the nominal problem without fixed cost, the optimal policy for the nomi- 
nal case is (S,S) with the threshold at time k being Sk = Wk- 

(c) For the nominal problem with fixed cost, if we denote by tj (j = 1, . . . , J) 
the times where stock is ordered and Sj , Sj the corresponding thresholds at time 
tj, we have: 



Ij 







Sj — r. w tj+ii 


(9) 






i — 0 




and 


ti-i 


Lj - 1 — 1 






Si = X 0 - ^2 Wi ’> 


Sj = - W L- i+i> 3 ^ 2 > 


(10) 




i = 0 


i—Ij — i+l 




where Lj 


II 

-bT 5 

+ 

ii 


pLj — c 1 {j=j} 
h+p 





We next present the main result regarding the structure of the optimal robust 
policy. 



Theorem 3 (Optimal robust policy). 

(a) The optimal policy in the robust formulation ( 8 ), evaluated at time 0 for 
the rest of the horizon, is the optimal policy for the nominal problem with the 
modified demand: 

w' k = Wk + P - k , ( A k - A k - 1 ) , (11) 

p + h 

where Ak = q^k + ]Ci=o r ik deviation of the cumulative demand from 

its mean at time k, q* and r* being the optimal q and r variables in (8). (By 
convention q - 1 = r .,-1 = 0.) In particular it is (S,S) if there is no fixed cost 
and (s, S) if there is a fixed cost. 
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(b) If there is no fixed cost, the optimal robust policy is (S,S) with S k = w' k for 
all k. 

(c) If there is a fixed cost, the corresponding thresholds Sj, Sj, where j = 
indexes the ordering times, are given by (9) and (10) applied to the modified 
demand w( 



k ’ 



(d) The optimal cost of the robust problem (8) is equal to the optimal cost for the 
nominal problem with the modified demand plus a term representing the extra 
cost incurred by the robust policy, 



%ph sr^T-l . 
L^k = 0 Ak ' 



n 



Proof. Let (u*, v*, q*, r*) be the optimal solution of (8). Obviously, setting the 
q and r variables to their optimal values q* and r* in (8) and resolving the linear 
programming problem will give u* and v* again. This enables us to focus on 
the optimal ordering policy only, taking the auxiliary variables q*, r* as given 
in the robust formulation (8). We have then to solve: 



T-l 

min E [cu k + Kl {uk>0} + max (h(x k+ i + A k ),p(-x k + 1 + H fe ))] (12) 

- k—0 



where x k +i = x 0 + J2i=o( u i ~ w i ) and A k = Q* k r k + o r ik for k - 

We define a modified stock variable x' k , which evolves according to the linear 
equation: 



x k+ 1 = x k + Uk - [ Wk + 



p — h 



(A k — A k -i) 



(13) 



with x' 0 = xo- Note that the modified demand w' k is not subject to uncertainty. 
We have: 



max (h(xk + 1 + A k ),p(-x k +i + A k )) = max (hx k+1 , 



-px' k+1 ) + 



2 ph 
p+h 



A k . (14) 



The reformulation of the robust model, given the optimal q* and r* variables, 
as a nominal inventory problem in the modified stock variable x k (plus the fixed 

cost o A k) follows from injecting (14) into (12). This proves (a) and 

(d). We conclude that (b) and (c) hold by invoking Lemma 1. □ 



Remark: For the case without fixed cost, and for the case with fixed cost when 
the optimal ordering times are given, the robust approach leads to the thresholds 
in closed form. For instance, if the demand is i.i.d. (w k = w, w k = w for all k), we 

have A k = w r k and, if there is no fixed cost, S k = w' k = w+ ^ . , w ( T k — A--i) 

P I Hj 

for all k. 

Hence, the robust approach protects against the uncertainty of the demand 
while maintaining striking similarities with the nominal problem, remains com- 
putationally tractable and is easy to understand intuitively. 
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3.2 The Capacitated Model 

So far, we have assumed that there was no upper bound either on the amount 
of stock that can be ordered or on the amount of stock that can be held in the 
facility. In this section, we consider the more realistic case where such bounds 
exist. The other assumptions remain the same as in Sect. 3.1. 

The Model with Capacitated Orders. The extension of the robust model 
to capacitated orders of maximal size d is immediate, by adding the constraint: 

Ufc < d, Vfc, (15) 

to (8) . We next study the structure of the optimal policy. 

Theorem 4 (Optimal robust policy). The optimal robust policy is the op- 
timal policy for the nominal problem with capacity d on the links and with the 
modified demand defined in (11). As a result, the optimal policy remains ( S,S ) 
(resp (s,S)) in the case without (resp with) fixed cost. 

The Model with Capacitated Inventory. We now consider the case where 
stock can only be stored up to an amount C. This adds the following constraint 
to (8): 

k 

Xo + ^2{ui - Wi) < C , (16) 

i= 0 

where Wi = Wi + w % • Zi such that z € {\zi\ < 1 Vi, Y)i = o \ z i\ — -He Vfc}. 
This constraint depends on the uncertain parameters Wi . Applying the technique 
developed in Sect. 2, we rewrite this constraint in the robust framework as: 

k 

x k + 1 + qk^k + ^2 r ik < C, Vfc, (17) 

i= 0 

where qu and r^ are defined in (8). We next study the optimal policy. 

Theorem 5 (Optimal robust policy). The optimal robust policy is the op- 
timal policy for the nominal problem subject to the modified demand defined in 
(11), and with inventory capacity at time 0 equal to C, and inventory capacity 
at time k + 1, k > 0, equal to C — A fc. 

As a residt, the optimal policy remains (S, S) (resp (s, S) ) in the case without 
(resp with) fixed purchasing cost. 
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4 The Network Case 

4.1 The Uncapacitated Model 

We now extend the results of Sect. 3 to the network case. We first study the case 
of tree networks, which are well suited to describe supply chains because of their 
hierarchical structure: the main storage hubs (the sources of the network) receive 
their supplies from outside manufacturing plants and send items throughout the 
network, each time bringing them closer to their final destination, until they 
reach the stores (the sinks of the network). Let S be the number of sink nodes. 
When there is only one sink node, the tree network is called a series system. 

We define echelon k, for k = 1, ... ,N with N the total number of nodes 
in the network, to be the union of all the installations, including k itself, that 
can receive stock from installation k, and the links between them. In the special 
case of series systems, we number the installations such that for k = 1 ,... ,N, 
the items transit from installation k + 1 to fc, with installation N receiving its 
supply from the plant and installation 1 being the only sink node, as in [7]. In 
that case, the demand at installation k + 1 at time t is the amount of stock 
ordered at installation k at the same time t. We also define, for k = 1, . . . , N: 
Ik{t) : the stock available at the beginning of period t at installation k, 

Xk ( t ) : the stock available at the beginning of period t at echelon k, 

D ik k(t) : the stock ordered at the beginning of period t at echelon k to its supplier 
iki 

W s (t) : the demand at sink node s during period t, s = 1, . . . , S. 

Let .ZV(fc) be the set of installations supplied by installation k and 0(k) the 
set of sink nodes in echelon k. We assume constant leadtimes equal to 0, backlog 
of excess demand, and linear dynamics for the stock at installation k at time 
t = 0,...,T-l : 



Ik(t + 1) — ifc(f) + Di k k(t) — Dkj(t), (18) 

jeN(k) 

By convention, if k is a sink node s, X^jeAf(fc) Dkj{t) = W s (t). This leads to the 
following dynamics for the stock at echelon k at time t = 0, . . . ,T — 1: 

X k (t+1) = X k (t) + D ihk (t) - W-W- ( 19 ) 

seO(fe) 

Furthermore, the stock ordered by echelon k at time t is subject to the 
coupling constraint: 



y Dki(t) < max(Jfc(f),0), Vfc, Vt, (20) 

ieJV(fc) 

that is, the total order made to a supplier cannot exceed what the supplier 
has currently in stock, or, equivalently, the supplier can only send through the 
network items that it really has. Since the network was empty when it started 
operating at time to = — oo, it follows by induction on t that Ik(t) > 0 for all 
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k > 2. Therefore the coupling constraint between echelons is linear and can be 
rewritten as: 

5] D ki (t)<X k (t)~ Y Vfc, Vt. (21) 

ieN(k) ieN(k) 

Finally, we specify the cost function. We assume that each echelon k has the 
same cost structure as the single installation modelled in Sect. 3.1 with spe- 
cific parameters ( Ck,K k ,h k ,Pk )• We also keep here the same notations and as- 
sumptions as in Sect. 3.1 regarding the uncertainty structure at each sink node. 
In particular, each sink node s has its own threshold sequence (t) evolving 
over time that represents the total budget of uncertainty allowed up to time 
t for sink s. We have W s (t) = W s (t) + W s {t) ■ Z s (t) such that the Z s (t ) be- 
long to the set V s = {\Z s (t)\ < 1 Vf, J2t=o ^s( T ) < r a (t), Vi}. We assume 
0 < r s (t) — r s (t — 1) < 1 for all s and t, that is, the budgets of uncertainty are 
increasing in t at each sink node, but cannot increase by more than 1 at each 
time period. 

Applying the robust approach developed in Sect. 2 to the holding/shortage 
constraints in the same manner as in Sect. 3, we obtain the mixed integer pro- 
gramming problem: 



min EE E {ckiDki{t) + KkiVki(t) + Yi(t)} 

t = o fe=iieJV(fc) 

s.t. Yi(t) > hi < Xi(t + 1) + Y + Y r s( T , t) J 

I sSO(i) V r=0 / 

Yi(t) >Pi < -Xi(t + 1) + Y [9s(t)r s (t) + Y r AY t) 
( seo(i) V t — o 

5] D ki (t) <X k (t) - Y 

ieJV(fe) ieAf(fe) 

q s (t) + r s (r, t) > W s (t), 
q s (t) > 0, r s (r, t) > 0, 

0 < D kl {t) < MV ki (t), V ki (t) G {0, 1}, 



with Xi(t- 1-1) = W(0) + Yt=o [ D ki(r) - J2 sgO(z) ^ / s( t )} for * and where 
k supplies i. 

As in the single-station case, an attractive feature of this approach is that 
the robust model of a supply chain remains of the same class as the nominal 
model, that is, a linear programming problem if there are no fixed costs and a 
mixed integer programming problem if fixed costs are present. Therefore, the 
proposed methodology is numerically tractable for very general topologies. The 
main result is as follows. 
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Theorem 6 (Optimal robust policy). 

(a) The optimal policy in (22) for echelon k is the optimal policy obtained, for 
the supply chain subject to the modified, deterministic demand at sink node s 
(for s e 0{k) ): 

W' S}k {t) = W s (t) + ^ (A s (t) - A s (t - 1)) , (23) 

Pk + h k 

where A s (t) = q*(t)T s (t) + r s ( T > ^)> an< ^ r s being the optimal q and r 

variables associated with sink node s in (22). 

(b) The optimal cost in the robust case for the tree network is equal to the optimal 
cost of the nominal problem for the modified demands, plus a term representing 
the extra cost incurred by the robust policy, 



E 



Zpkhk 
Pk + hk 



T - 1 



E E 

t—0 sGO(fe) 



The case of more general supply chains is complex because they cannot be 
reduced to a tree network: the need might arise to order from a more expensive 
supplier when the cheapest one does not have enough inventory. We can still de- 
fine echelons for those networks in a similar manner as before, and the evolution 
of the stock at echelon k, which is supplied by the set of installations I ( k ) and 
has the set 0{k) as its sink nodes, is described by the following linear equation: 



X k {t + 1) = X k (0) + E I E - E W i( T ) \ ■ ( 24 ) 

r=o [iei(fc) j'eo(fc) J 

With the standard cost assumptions used before, the echelons cannot be studied 
independently and the optimal policy is not necessarily basestock, even in the 
simple case of demand without uncertainty. This is illustrated by the following 
example. 




Fig. 1. A network for which the optimal policy is not basestock. 



The network in Fig. 1 has two installations, and therefore two echelons. 
Echelon 1 can be supplied by installation 2 at a unit cost ci = 1, without 
any fixed ordering cost, and has the option to order directly from the plant for 
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the same unit cost C 2 = 1, but with an additional fixed cost K 2 = 4 incurred 
whenever an order is made. This option is attractive only if installation 2 does 
not have enough stock in inventory. The holding and shortage unit costs at 
echelon 1 are h\ = p\ = 2. The horizon is 1 time period, and the demand at 
time 0 is deterministic, equal to 10 units. Echelon 1 has only 5 units in inventory 
at time 0. 

Comparing the two options, it is easy to see that it is optimal for echelon 1 
to order 5 units from 2 if 2 has 5 units in inventory at time 0, 2 units from 2 
and none from the plant if 2 has 2 units, and 5 units from the plant if 2 has no 
item in inventory. Therefore, the optimal amount of stock on hand and on order 
at echelon 1 at time 0 is 10, resp. 7, 10, units if installation 2 has 5, resp 2, 0 
units in inventory at time 0. Thus, the optimal policy is not basestock. 

Also, while we can reformulate the robust problem as a new problem with 
modified demand in the same fashion as before, it loses some of its meaning since 
distinct echelons can now “see” the same sink node but different demands at this 
node (because of the cost parameters specific to each echelon, which appear in 
the expression of the modified demand). Hence, it is not clear how they can work 
together to meet this demand optimally. 

However, the proposed robust methodology remains numerically tractable 
in a wide range of settings, in particular with a holding/shortage cost at the 
installation level instead of the echelon. This illustrates the applicability of the 
proposed approach to different cost structures. 

4.2 The Capacitated Model 

We now refine our description of the inventory problem in a supply chain by 
introducing upper bounds on the amount of stock that can be ordered and/or 
held at any time and at any echelon. As explained in Sect. 3.2, an upper bound 
on the maximal order can be directly introduced in the proposed approach, by 
adding the constraint: 

D ki (t) < d ki Vfc, Vi e Af(k% Vi, (25) 

to (22). Inventory capacity, however, requires further manipulation, since the 
level of inventory held at an echelon at any time depends on the demand, which is 
subject to uncertainty. Similar manipulations as in Sect. 3.2 lead to the constraint 
Vfc, Vi: 

X k (t + 1)+ ^2 ( 9s(i)^ s (i) + ^r s (r,i) J < (26) 

seO(k) V r=0 / 

to be added to the formulation, q(t) and r(r, i) being defined as in (23). 

We next study the structure of the optimal policy. 

Theorem 7. The optimal policy at each echelon remains basestock in presence 
of link and echelon capacities. It is identical to the optimal policy of a nominal 
problem at a single station subject to the modified demand defined in (23), time- 

varying echelon capacities: C' k (t + 1) = C k — +h k ^-'seO(fc) A s {t), where C k is 
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the original capacity at echelon k, and link capacities that incorporate Dki(t) < 
dki for all k, i £ Af(k) and t, as well as the capacity induced by the coupling 
constraint (21). 



5 Numerical Implementation 



We now apply the proposed methodology to the example of minimizing the cost 
at a single station. The horizon is T = 10 time periods, the initial inventory is 
Xo = 150, with an ordering cost per unit c = 1, a holding cost h = 2 and a short- 
age cost p = 3, in appropriate measurement units. There is no fixed ordering cost. 
The stochastic demand is i.i.d. with mean w = 100. In the robust framework, 
we consider that the demand belongs to the interval [0,200], that is w = 100. 
We compare the expected costs of the robust policy and of the stochastic policy 
obtained by dynamic programming as a function of the standard deviation a of 
the distribution. 

To select the budgets of uncertainty we minimize a tight upper bound on the 
expected cost of the system over the set of nonnegative distributions of given 
mean and variance, using the results in [4]. This yields for all k: 



r/. = min 







k + 1 




and the modified demand at time k is in this example: 

w' k = w+ ° a (Vk + 1 - Vk), 
V 1 — or 



(27) 



(28) 



with a 



' 

p + h' 



Expected costs are computed using the mean of a sample of 



size 1,000. 

In the first set of experiments, the stochastic policy is computed using a 
binomial distribution. In the second set of experiments, the stochastic policy is 
computed using an approximation of the gaussian distribution on seven points 
(w — 3cr, w — 2cr, . . . , w + 2cr, w + 3(7 ) . In both cases, the actual distribution is 
Gamma, Lognormal or Gaussian, with the same mean W and standard deviation 
cr. The impact of the mistake on the demand distribution is measured by the 
ratio (DP — ROB) / DP , with DP , resp. ROB , the expected cost obtained using 
dynamic programming, resp. the robust approach. The results are shown in Fig. 
2. In the first case, where the distributions are very different beyond their first 
moments, the impact of the ratio increases as the standard deviation increases 
and the robust policy outperforms dynamic programming by up to 8%. In the 
second case, the two methods are equivalent in terms of performance, since the 
robust policy outperforms dynamic programming by at most 0.3%, which is not 
statistically significant. 

In Figs. 3-5, we study the impact of the cost parameters c, h and p in the 
settings described above, where we vary one parameter at a time. The impact 
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Fig. 2. Impact of the standard deviation on performance. 



of a change in the parameters is qualitatively similar in both cases, with lit- 
tle dependence on the actual distribution of the demand. The robust approach 
outperforms the stochastic policy for a wide range of parameters, although the 
stochastic policy leads to better results for large values of the ratio p/h (greater 
than about 3). The exact numbers depend on the distribution used to compute 
the stochastic policy. 




s 

i 



— Gamma 

- - Lognormal 



Fig. 3. Impact of the ordering cost. 



Overall the numerical evidence suggests that the robust policy performs sig- 
nificantly better than dynamic programming when assumed and actual distri- 
butions differ widely despite having the same mean and standard deviation, and 
performs similarly to dynamic programming when assumed and actual distribu- 
tions are close. The results are thus quite promising. 
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Fig. 4. Impact of the holding cost. 





Fig. 5. Impact of the shortage cost. 



6 Conclusions 

We have proposed a deterministic, numerically tractable methodology to ad- 
dress the problem of optimally controlling supply chains subject to uncertain 
demand. Using robust optimization ideas, we have built an equivalent model 
without uncertainty of the same class as the nominal problem, with a modi- 
fied demand sequence. Specifically, the proposed model is a linear programming 
problem if there are no fixed costs throughout the supply chain and a mixed 
integer programming problem if fixed costs are present. 

The key attractive features of the proposed approach are: (a) It incorpo- 
rates a wide variety of phenomena, including demands that are not identically 
distributed over time and capacity on the echelons and links; (b) it uses very 
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little information on the demand distributions; (c) it leads to qualitatively sim- 
ilar optimal policies (basestock policies) as in dynamic programming; (d) it is 
numerically tractable for large scale supply chain problems even in networks, 
where dynamic programming methods face serious dimensionality problems; (e) 
in preliminary computational experiments, it often outperforms dynamic pro- 
gramming based solutions for a wide range of parameters. 



References 

1. Ben-Tal, A., Nemirovski, A. (2000): Robust solutions of linear programming prob- 
lems contaminated with uncertain data, Math. Programm. , Ser. A 88:411-424- 

2. Ben-Tal, A., Nemirovski, A. (1999): Robust solutions of uncertain linear programs, 
Oper. Research Lett. 25, 1-13. 

3. Ben-Tal, A., Nemirovski, A. (1998): Robust convex optimization, Math. Oper. Res. 
23, 769-805. 

4. Bertsimas, D., Popescu, I. (2002): On the relation between option and stock prices: 
a convex optimization approach, Oper. Res. 50. 358-374 ■ 

5. Bertsimas, D., Sim, M. (2003): The price of robustness, to appear in Operations 
Research. 

6. Bertsimas, D., Sim, M. (2003): Robust discrete optimization and network flows, 
Mathematical Programming , Vol. 98, No. 1-3, 49-71. 

7. Clark, A., Scarf, H. (1960): Optimal policies for a multi-echelon inventory problem, 
Management Science Vol. 6 No. 4, 475-490. 

8. El-Ghaoui, L., Lebret, H. (1997): Robust solutions to least-square problems to 
uncertain data matrices, SIAM J. Matrix Anal. Appl. 18, 1035-1064 ■ 

9. El-Ghaoui, L., Oustry, F., Lebret, H. (1998): Robust solutions to uncertain semidef- 
inite programs, SIAM J. Optim. 9, 33-52. 

10. Scarf, H. (1958): A min-max solution of an inventory problem, Studies in the 
mathematical theory of inventory and production, 201-209. 

11. Soyster, A.L. (1973): Convex programming with set-inclusive constraints and ap- 
plications to inexact linear programming, Oper. Res. 21. 1154-1157. 

12. Zipkin, P. (2000): Foundations of inventory management, McGraw-Hill Higher Ed- 
ucation. 




Hedging Uncertainty: Approximation 
Algorithms for Stochastic Optimization 
Problems* 



R. Ravi and Amitabh Sinha 

Graduate School of Industrial Administration, Carnegie Mellon University, 
Pittsburgh, PA, USA 
ravi@cmu . edu, asinha@andrew . emu . edu 



Abstract. We study the design of approximation algorithms for stoch- 
astic combinatorial optimization problems. We formulate the problems in 
the framework of two-stage stochastic optimization, and provide nearly 
tight approximations. Our problems range from the simple (shortest 
path, vertex cover, bin packing) to complex (facility location, set cover), 
and contain representatives with different approximation ratios. 

The approximation ratio of the stochastic variant of a typical problem is 
of the same order of magnitude as its deterministic counterpart. Further- 
more, common techniques for designing approximation algorithms such 
as LP rounding, the primal-dual method, and the greedy algorithm, can 
be carefully adapted to obtain these results. 



1 Introduction 

With the increasing success of optimization algorithms in process optimization, 
these methods are making inroads into earlier planning stages of large scale 
projects. The inherent difference between optimization at the planning stage 
and post-facto optimization is that in the former, data is not fully available. 
Yet costly decisions with wide-ranging implications need to be taken in the face 
of incomplete data. Nevertheless, quite often forecasts of future uncertainty are 
available that can be used in the planning model. Forecasts, by nature, are 
imprecise and provide at best a range of possible futures. The field of stochastic 
optimization is an attempt to model such situations. For a detailed introduction, 
the reader is referred to one of the recent texts on the topic [4,19]. 

In a parallel development, the field of approximation algorithms [33,2] evol- 
ved to counter the prohibitive resource requirements for exact solution of NP- 
lrard combinatorial optimization problems. Informally, these algorithms run in 
polynomial time and deliver a performance ratio on the quality of the output 
solution over all instances. As the size of the models being solved increases in 
scale, this solution approach gains in importance. 

* This work was supported in part by NSF grant CCR-0105548 and ITR grant CCR- 
0122581 (The ALADDIN project). 
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However, as approximation algorithms become more sophisticated in scope 
and technique, the refrain from real-world practitioners who have the need for 
such algorithms is that the input data is seldom well-defined, thus diminishing 
the value of the solutions and guarantees provided by the algorithm. Conversely, 
while the Held of stochastic optimization models the uncertainty in data fairly 
well, the running times of the exact algorithms developed in the stochastic op- 
timization community often prove prohibitive. This paper combines the best of 
both worlds, by providing approximation algorithms for the stochastic version 
of several classical optimization problems. 

2 Background 

Two- stage Model Among the most popular models in stochastic optimiza- 
tion is the two-stage model with recourse. At the outset, some data may be 
known deterministically, whereas the uncertain future is characterized only by a 
probability distribution. The decisions made at this point are referred to as the 
first-stage decisions. Subsequently, the actual future is realized, and then there 
may be the opportunity to augment the first-stage solution in order to optimize 
for the realized scenario. This second stage of decision making is called the re- 
course stage. The goal is to optimize the first-stage decision variables so as to 
minimize the expected cost over both stages. 

Mathematical Formulation We consider the two-stage stochastic optimiza- 
tion problem with recourse, with an additional restriction: finitely many sce- 
narios. This means that the future will be one of a finite set of possibilities 
(scenarios) , and the parameters and probability of occurence of each scenario is 
known up-front. The mathematical formulation of this model is given below. 

Vector x° is the set of first-stage decision variables, with constraints Ax° — b , 
and a cost vector is c. There are m scenarios, with the k th scenario having 
probability of occurence pk, cost vector q k , and decision variables x k . If the 
k th scenario is realized, then the combined solution ( x°,x k ) must satisfy the 
constraints given by the matrix T k and requirement vector h k . Let P denote 
additional constraints such as non-negativity or integrality. 

min c T x° + Pk(q k ) T x k (IPs) 

s.t. Ax° = b 

T k (x°, x k ) — h k k = 1, 2, . . . , m 

(a: 0 , x k ) € P k = 1, 2, . . . , m 

The interested reader may refer to any of the texts cited above for a more 
complete description of models of stochastic optimization and their uses. Schultz, 
Stougie and van der Vlerk [30] provide an excellent survey of two-stage stochastic 
integer programming, while Kong and Schaefer [20] recently provided approxi- 
mation algorithms for a class of such problems. 

Approximation Algorithms The raison d’etre for approximation algorithms 
is the prohibitively high running time of exact algorithms for integer program- 
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Problem 


Det. 

approx. 


Stochastic 

elements 


Our results 
(Apx. ratio) 


Hardness 


Bin packing 


APTAS [5] 


Object sizes 


APTAS 


NP-comp. [5] 


Shortest paths 


1[7] 


Sink only 
Sink and metric 


5 

0(log 2 nlog m) 


MAX-SNP 
l?(log 2 n) 


Vertex cover 


2 [27] 


Vertex weights, incidence 


2 


1.16[14] 


Facility location 


1.52 [25] 


Demands, facility costs 


8 


1.46 [10] 


Set cover 


0(logn)[17] 


Set weights, inclusions 


0(log nm) 


!7(logn)[l] 
12(log m) 



Fig. 1 . Summary of results. We use m to denote the number of scenarios and n 
to refer to the number of combinatorial elements (number of vertices in graph 
problems and number of elements in the set cover problem) . 



ming and combinatorial optimization, due to their NP-completeness. Approx- 
imation algorithms run in time polynomial in the size of the input, and also 
provide guarantees in the form of approximation ratios. If C is a class of min- 
imization problems and OPT {I) and A(I) respectively denote the value of an 
optimal solution and the solution output by algorithm A, then the approximation 
ratio p(A) of algorithm A is defined as p(A) = max/gc opt } i ) • 

While the area of approximation algorithms has been a very active field, most 
approximation algorithms assume complete knowledge of the input at the out- 
set, barring a few exceptions such as scheduling problems [26,32]. Recently, and 
independently of us, Immorlica et al. [15] considered approximation algorithms 
for stochastic optimization problems with a restriction on the cost function: in 
the second stage, all costs go up uniformly by a factor of A. A generalization of 
their model was considered by Gupta et al. [12], who provided approximation 
algorithms with only sampling access to the second-stage realization process, 
thereby obtaining a framework which can handle an arbitrary number of scenar- 
ios as long as there is an efficient process to sample the second stage. 

Our Results We demonstrate the relevance and applicability of developing ap- 
proximation algorithms for stochastic optimization. We carefully adapt existing 
techniques for deterministic versions of several problems to provide approxima- 
tion guarantees for the stochastic versions within constant factors. 

Our results are summarized in Figure 1. The current best known deterministic 
approximations are listed, with a “1” meaning that the problem can be solved 
optimally in polynomial time. All the stochastic approximation ratios are derived 
in this paper. Some of the hardness results are carried over from the underlying 
deterministic problems; the remaining are proved in this paper. An APTAS is an 
asymptotic polynomial time approximation scheme, which is an algorithm whose 
performance ratio approaches 1 as the number of objects increases. A problem 
is said to be MAX-SNP-lrard [28], abbreviated MAX-SNP in the table, if there 
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is a constant c > 1 such that it is impossible to approximate the problem with 
performance ratio smaller than c unless P = NP. 

Paper Outline In the sequel, we consider the five problems listed in Figure 1. 
We consider Stochastic Vertex Cover in the next section, and prove our results. 
Subsequent sections are devoted to the stochastic versions of four other problems 
that we study: Facility Location, Shortest Paths, Bin Packing and Set Cover. 
We conclude with directions for future research in Section 8. 

3 Vertex Cover 

We are given a first-stage (undirected) graph G = (V,E 0 ), with m possible 
scenarios, each consisting of a probability of occurrence pk and a set of edges Ek 
(not necessarily subsets of E 0 ). The first-stage cost of vertex v is c°, and its cost 
in scenario k is c^. The objective is to identify a set of vertices to be selected in 
the first stage, so that the expected cost of extending this set to a vertex cover 
of the edges of the realized second-stage scenario is minimized. 

We require that the edges in Ek \ E 0 have to be covered in the second stage, 
even if one of their end-points was chosen in the first stage. This generalizes 
the case when first-stage vertices cover all second-stage edges incident to them. 
The best known approximation algorithm for the deterministic version of vertex 
cover has performance ratio 2 — ^log^jq ? due to Monien and Speckenmeyer 
[27]. A lower bound of 1.16 on the hardness of approximating the problem was 
shown by Hastad [14]. Our algorithm for the generalized stochastic version of 
vertex cover has approximation ratio 2, asymptotically matching the best known 
approximation for the deterministic version. 

Integer Program Formulation In formulation IPsvc below, variable x * 
indicates whether or not vertex v is purchased in scenario k (where k = 0 denotes 
the first stage) . Edges in Ek D Eq may be covered in either the first or second 
stage, while edges in Ek \ E 0 must be covered in the second stage. Observe that 
a 4-approximation can be easily obtained by rounding IPsvc [15]- We obtain a 
tighter approximation ratio using a primal-dual algorithm. 

\E\\\c°x° + YJk=\PkC k x k (IPsvc) 
s.t. x ° + x° + x’l + > 1 Vuv G E k fl E 0 ,\/k 

x 1 ^ + x^, > 1 Muv G E k \ E 0 yk 

x non-neg. integers 

Dual Program The dual of the linear relaxation of IPsvc is shown below. 
Variable packs edge e in Ek if e G Ek , and it packs e G Eo if e G Ek fl E$. 

max X)JT=i J2u,vev Vuv ( DP SV c ) 

S-t. T,eGE k :vee Ve < Pk 4 V«, Vfc 
X)fc=l ^2eeE 0 nE k :vee Ve ^ c v ^ v 

y> o 
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Algorithm The algorithm is a greedy dual-ascent type of primal-dual algo- 
rithm, with two phases. In Phase I, we raise the dual variables y k uniformly 
for all edges in E k \ Eq, separately for each k from 1 to m. All vertices which 
become tight (have the first dual constraint packed to p k c k ) have x k set to 1, 
and are deleted along with adjacent edges. We proceed this way until all edges 
in E k \ E° are covered and deleted. 

In Phase II, we do a greedy dual-ascent on all uncovered edges of E k , which 
are contained in E k r\E 0 . We raise y k for all uncovered edges for k = 0 to m. We 
use a slightly different rule for purchasing vertices: If a vertex is tight for x® (i.e. , 
second dual constraint packed to c°), then we select it in the stage 1 solution by 
setting = 1, and if it is not tight for x° but is tight for x k (packed in the first 
dual constraint), then we select it in the recourse solution and set x k = 1. 

Theorem 1. The integer program IPsvc can be rounded by the primal-dual 
algorithm described above within a factor of 2 in polynomial time. 

Proof. Consider an edge e = uv in scenario k. We must have selected one of its 
two end-points in either Phase I or Phase II (or both), so the algorithm yields a 
feasible solution. We use linear programming duality to bound the cost of the so- 
lution by showing that the cost of our solution is no more than 2 „ e y Vuvi 

where y is the dual solution constructed by our algorithm. Each time we set an 
x k variable to 1, we assign some dual variables to it such that (i) the sum of 
dual variables assigned to each such x k variable equals puc k (where po = 1), and 
(ii) each dual variable is assigned at most twice. 

Consider a vertex v which was selected (ie, x k was set to 1) in scenario k 
in either Phase I or Phase II. We assign all dual variables y k such that v is an 
end point of e to this vertex, and since v is selected only when the constraint 
J2eeE k :v£e Ve < PkC k goes tight, we maintain (i). An edge e in E k \E 0 is assigned 
to a vertex v only if x k is set to 1 for k ^ 0, and since an edge has at most 2 
end-points, we ensure (ii) for edges in E k . 

Next consider a vertex v for which we set x'l to 1. Therefore, the constraint 
ElEe &Eo nE k :veeye < c° must have gone tight, and all edges in the sum are 
assigned to the variable x® . This assignment once again ensures (i) . This assign- 
ment only includes edges in E 0 C\ E k , and these edges are not assigned to any 
variable x k for k ^ 0, thus also ensuring (ii) for all edges in E 0 fl E k . 

These two cases cover all the possibilities, thus proving the theorem. 

4 Facility Location 

As in the classical uncapacitated facility location problem, we are given a set 
of facilities F and a set of clients D, with a metric Cjj specifying the distances 
between every client and every facility. However, the demand of each client is 
not known at the first stage. In scenario k, client j has demand d k , which may 
be zero. Facility i has a first-stage opening cost of ff, and recourse costs of f k 
in scenario k. These may be infinity, reflecting the unavailability of the facilities 
in various scenarios. We abbreviate this problem as SFL. 
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m 

min Y fiVi + Y Pk ( fi y i + d hA) ( IP SFL ) 

iFF fc= 1 ieF iGFj'S-D 

S.t. Yl^o - Vj G D,V/c 

< Zh 9 + Vi v« G F, Vj G D, Vfc 

x, y non-negative integers 

The problem is best explained by the IP formulation IPsfl above. While 
our algorithm extends to arbitrary demands, for simplicity we only study the 
case when all d k, s are either 0 or 1. Variable is 1 if and only if client j is 
served by facility i in scenario k. If x p = 1, then facility i must either be opened 
at the first stage ( y 9 = 1) or in recourse in scenario k (y k = 1) (or both). 

History and Non-triviality of the Problem The classical (determinis- 
tic) uncapacitated facility location problem has a rich history (see Cornuejols, 
Nemhauser and Wolsey [6] for a survey). Balinski [3] introduced an integer pro- 
gramming formulation for this problem which has led to several approximation 
algorithms. The first constant factor approximation using this formulation is 
due to Shmoys, Tardos and Aardal [31], and the current best algorithm, due 
to Mahdian, Ye and Zhang [25] uses a formulation which differs only slightly. 
Indeed, our formulation (IPsfl) extends Balinski’s formulation to the stochas- 
tic setting. In the stochastic optimization community, Louveaux and Peeters [23] 
considered a slightly different version of stochastic facility location, and provided 
a dual-ascent based exact (non-polynomial time) algorithm for it. 

Notice that if the second stage facility costs were identical to those in the first 
stage for all scenarios, then we can “de-couple” the stochastic components of the 
problem and solve for each scenario independently. On the other hand, if there 
was no second stage and all facilities had to be opened in the first stage, then 
SFL reduces to an instance of the usual UFL, where the probability multipliers in 
the expected service costs can be incorporated into the demand terms (thinking 
of the demands as being scaled based on the probability of occurrence); in this 
case, existing approximations for UFL apply directly. 

The added difficulty, and indeed the interesting aspect of the model, arises 
from varying (and typically increased) second-stage facility costs under different 
scenarios. We cannot treat each scenario by itself, since the different scenarios 
interact in utilizing first-stage facilities. Our algorithm has to carefully balance 
the effects of the two stages in deciding what facilities to open. 

Algorithm Our approximation algorithm proceeds along the lines of the LP- 
rounding algorithm of Shmoys, Tardos and Aardal [31], with some crucial differ- 
ences. We begin by solving the linear relaxation of IPsfl ■ Let ( x,y ) denote an 
optimal LP solution. The first step in rounding this fractional solution is using 
the filtering technique of Lin and Vitter [22]. We fix a constant 0 < a < 1. For 
every client-scenario pair (j, fc), we define its optimal fractional service cost to 
be c)j k = Y^i°ij x ij- Order the facilities which serve the pair (j, k) according to 
non-decreasing distance from j. The a point 5 yfc(ct) for the client-scenario pair 




